
#import basic dependencies (import1)
import pandas as pd
import numpy as np

import datetime
import math

import matplotlib.pyplot as plt
from matplotlib import rcParams
import seaborn as sns

from sklearn.preprocessing import StandardScaler

import warnings
#warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.filterwarnings("ignore")


df_full = pd.read_csv("../data/external/feature_data_substation_bancroft_labelled.csv")
df_full.head()


print(df_full.info())
print(df_full.dtypes)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 65345 entries, 0 to 65344
Columns: 184 entries, Unnamed: 0 to month
dtypes: float64(166), int64(8), object(10)
memory usage: 91.7+ MB
None
Unnamed: 0          int64
substation         object
run_datetime       object
valid_datetime     object
horizon           float64
                   ...   
label3              int64
wep                object
storm_id2          object
year                int64
month               int64
Length: 184, dtype: object


# Change datetime columns to datetime format
df_full['run_datetime'] = pd.to_datetime(df_full['run_datetime'])
df_full['valid_datetime'] = pd.to_datetime(df_full['valid_datetime'])

# Add temperature column for degree in celsius
df_full['avg_temp_celsius'] = df_full['avg_temp'] - 273.15

# Translate wind directions into cardinal directions
def wind_direction_label(degrees):
    if 22.5 < degrees <= 67.5:
        return 'Northeast'
    elif 67.5 < degrees <= 112.5:
        return 'East'
    elif 112.5 < degrees <= 157.5:
        return 'Southeast'
    elif 157.5 < degrees <= 202.5:
        return 'South'
    elif 202.5 < degrees <= 247.5:
        return 'Southwest'
    elif 247.5 < degrees <= 292.5:
        return 'West'
    elif 292.5 < degrees <= 337.5:
        return 'Northwest'
    else:
        return 'North'
    
df_full['wind_direction_label'] = df_full['avg_winddir'].apply(wind_direction_label)


data = {
    'Name': [
        'run_datetime',
        'wep',
        'avg_temp',
        'min_wet_bulb_temp',
        'avg_dewpoint',
        'avg_temp_change',
        'avg_windspd',
        'max_windgust',
        'avg_winddir',
        'wind_direction_label',
        'max_cumulative_precip',
        'max_snow_density_6',
        'max_cumulative_snow',
        'max_cumulative_ice',
        'avg_pressure_change'
        ],
    'Description': [
        'Date and time when the weather observations were recorded.',
        'Weather Event Type (WEP) is a categorization of weather conditions based on specific parameters at a given location and time.',
        'The average temperature measured at two meters above ground level, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Minimum wet bulb temperature recorded during the observation period, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Average dewpoint temperature observed during the recording period, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Average change in temperature during the observation period, obtained by calculating the difference between this observation and the following.',
        'Average wind speed measured during the recording period, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Maximum wind gust observed during the recording period, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Average wind direction (in degree) observed during the recording period, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Wind direction (in cardinal direction) observed during the recording period, obtained by recalculating the avg_winddir parameter.',
        'Maximum cumulative precipitation recorded, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Maximum snow density at a depth of 6 inches, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Maximum cumulative snow recorded, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Maximum cumulative ice recorded, considering all sensors, for the entire duration of one hour in the Bancroft region.',
        'Average change in atmospheric pressure during the observation period, considering all sensors, for the entire duration of one hour in the Bancroft region.'
        ],
    'Role': [
        'ID / predictor', 
        'response', 
        'response / predictor', 
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor',
        'predictor'
        ],
    'Type': [
        'numerical continuous / ID', 
        'categorical nominal', 
        'numerical continuous', 
        'numerical continuous', 
        'numerical continuous', 
        'numerical continuous', 
        'numerical continuous', 
        'numerical continuous', 
        'numerical continuous', 
        'categorical ordinal', 
        'numerical continuous',
        'numerical continuous',
        'numerical continuous',
        'numerical continuous',
        'numerical continuous'
        ],
    'Format': [
        type(df_full['run_datetime'][0]),
        type(df_full['wep'][0]), 
        type(df_full['min_wet_bulb_temp'][0]),
        type(df_full['avg_temp'][0]), 
        type(df_full['avg_dewpoint'][0]), 
        type(df_full['avg_temp_change'][0]), 
        type(df_full['avg_windspd'][0]), 
        type(df_full['max_windgust'][0]),
        type(df_full['avg_winddir'][0]), 
        type(df_full['wind_direction_label'][0]), 
        type(df_full['max_cumulative_precip'][0]), 
        type(df_full['max_snow_density_6'][0]),
        type(df_full['max_cumulative_snow'][0]), 
        type(df_full['max_cumulative_ice'][0]), 
        type(df_full['avg_pressure_change'][0])
]
}

data_dict_df = pd.DataFrame(data)

# Display the data dictionary
data_dict_df


#Creating a new dataframe with only the used columns (saves resources)
columns_to_keep = [
    'run_datetime',
    'wep',
    'avg_temp',
    'avg_temp_celsius',
    'min_wet_bulb_temp',
    'avg_dewpoint',
    'avg_temp_change',
    'avg_windspd',
    'max_windgust',
    'avg_winddir',
    'avg_winddir_sin',
    'avg_winddir_cos',
    'wind_direction_label',
    'max_cumulative_precip',
    'max_snow_density_6',
    'max_cumulative_snow',
    'max_cumulative_ice',
    'avg_pressure_change',
    'label0', #additional wep labels for various blue sky events
    'label1', #additional wep labels to differentiate between blue sky and extreme weather events
    'label2'
]

# Use loc to select only the specified columns
df = df_full.loc[:, columns_to_keep]
df


from sklearn.model_selection import StratifiedShuffleSplit

# Initialize StratifiedShuffleSplit
sss = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=42)

# Split the data using StratifiedShuffleSplit
for train_index, test_temp_index in sss.split(df, df['wep']):
    train_data = df.iloc[train_index]
    temp_data = df.iloc[test_temp_index]

# Further split temp_data into test and validation sets
for test_index, validation_index in sss.split(temp_data, temp_data['wep']):
    test_data = temp_data.iloc[test_index]
    validation_data = temp_data.iloc[validation_index]

# Display the shapes of the resulting sets
print("Train Data Shape:", train_data.shape)
print("Test Data Shape:", test_data.shape)
print("Validation Data Shape:", validation_data.shape)

# Display the distribution of labels 'wep'
print("\nDistribution of 'wep' in Train Data:")
print(train_data['wep'].value_counts())

print("\nDistribution of 'wep' in Test Data:")
print(test_data['wep'].value_counts())

print("\nDistribution of 'wep' in Validation Data:")
print(validation_data['wep'].value_counts())

Train Data Shape: (52276, 21)
Test Data Shape: (10455, 21)
Validation Data Shape: (2614, 21)

Distribution of 'wep' in Train Data:
wep
Blue sky day                                           42106
Mild snowfall                                           3598
Moderate snowfall                                       2336
Moderate rain                                           2104
Heavy snowfall with accumulated snow                     832
Frondurchlauf / Continuous freezing rain                 747
Storm with freezing rain / Heavy snow- and icestorm      344
Snowstorm with high precipitation                        209
Name: count, dtype: int64

Distribution of 'wep' in Test Data:
wep
Blue sky day                                           8421
Mild snowfall                                           719
Moderate snowfall                                       467
Moderate rain                                           421
Heavy snowfall with accumulated snow                    166
Frondurchlauf / Continuous freezing rain                150
Storm with freezing rain / Heavy snow- and icestorm      69
Snowstorm with high precipitation                        42
Name: count, dtype: int64

Distribution of 'wep' in Validation Data:
wep
Blue sky day                                           2106
Mild snowfall                                           180
Moderate snowfall                                       117
Moderate rain                                           105
Heavy snowfall with accumulated snow                     42
Frondurchlauf / Continuous freezing rain                 37
Storm with freezing rain / Heavy snow- and icestorm      17
Snowstorm with high precipitation                        10
Name: count, dtype: int64


df.describe().T


import matplotlib.dates as mdates

rcParams['figure.figsize'] = 20, 15
sns.set(style="darkgrid")

parameters = [
        ('avg_temp', 'Kelvin'),
        ('min_wet_bulb_temp', 'Kelvin'),
        ('avg_dewpoint', 'Kelvin'),
        ('avg_temp_change', 'Kelvin'),
        ('avg_windspd', 'm/s'),
        ('max_windgust', 'm/s'),
        ('avg_winddir', 'degrees'),
        ('max_cumulative_precip', 'mm'),
        ('max_snow_density_6', 'kg/m^3'),
        ('max_cumulative_snow', 'mm'),
        ('max_cumulative_ice', 'mm'),
        ('avg_pressure_change', 'Pa')
    ]

def calculate_moving_averages(df, variables, time_periods): # has already been done in previous work at IBM (by Julian Erath) for the dataset # selected_variables = ['max_windgust', 'min_wet_bulb_temp', 'max_snow_density_6', 'max_cumulative_precip', 'max_cumulative_ice', 'max_cumulative_snow', 'avg_windspd','avg_pressure_change','avg_temp','avg_dewpoint','avg_winddir','avg_temp_change'] # selected_time_periods_week = [168] # selected_time_periods_month = [720]
    for variable in variables:
        for period in time_periods:
            ma_column_name = f'MA_{variable}_{period}_hours'
            df[ma_column_name] = df[variable].rolling(period, min_periods=1).mean()



def plot_weather_parameters(ax, data, parameter, unit):
    ax.plot(data.run_datetime, data[parameter], label=f'{parameter} Actual')
    ax.plot(data.run_datetime, data[f"MA_{parameter}_week"], label=f'{parameter} Weekly MA')
    ax.plot(data.run_datetime, data[f"MA_{parameter}_month"], label=f'{parameter} Monthly MA')
    ax.set_title(f'Bancroft {parameter}')
    ax.set_ylabel(unit)

    # Set the x-axis major locator to YearLocator and the formatter to the year
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

    # Manually add the mid-year tick for July 1st
    mid_years = [mdates.date2num(datetime.datetime(year, 7, 1)) for year in range(data.run_datetime.dt.year.min(), data.run_datetime.dt.year.max() + 1)]
    ax.set_xticks(list(ax.get_xticks()) + mid_years, minor=True)
    ax.xaxis.set_minor_formatter(mdates.DateFormatter('/%m'))

    # Rotate the labels for better visibility
    plt.setp(ax.xaxis.get_majorticklabels(), rotation=90, ha='right')

def create_subplots(data):
    fig, axes = plt.subplots(4, 3, figsize=[20,15])

    for i, (parameter, unit) in enumerate(parameters):
        plot_weather_parameters(axes[i // 3, i % 3], data, parameter, unit)

    plt.suptitle('Bancroft weather parameters time series', fontsize=24)
    plt.tight_layout()
    # Update the legend to handle multiple lines
    handles, labels = axes[0, 0].get_legend_handles_labels()
    new_labels = [label.replace("avg_temp", "Average Temp.").replace("MA", "Moving Average") for label in labels]
    fig.legend(handles, new_labels, loc='upper right', ncol=3)

    plt.show()
    plt.savefig('../output/Bancroft weather parameters time series.png')

create_subplots(df_full)

<Figure size 2000x1500 with 0 Axes>


def plot_histogram(data, parameter, unit, ax):
    # Plot the histogram
    sns.histplot(data[parameter], bins=30, kde=True, ax=ax, color='blue')
    # Calculate the mean of the parameter
    median_value = data[parameter].median()
    # Draw a vertical line for the median value
    ax.axvline(median_value, color='red', linestyle='dashed', linewidth=2)
    # Add a label for the line
    ax.text(median_value, ax.get_ylim()[1]*0.9, f'Median: {median_value:.2f}', color='red')

    # Set title and axis labels
    ax.set_title(f'Histogram for {parameter}')
    ax.set_ylabel('Frequency')
    ax.set_xlabel(unit)

# Set up the subplots with shared y-axis
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=[20,15])

# Plot histograms for all parameters
for i, (parameter, unit) in enumerate(parameters):
    row, col = divmod(i, 3)
    plot_histogram(df, parameter, unit, axes[row, col])

# Manually set bottom margin to make x-axis visible
plt.subplots_adjust(bottom=0.1)  # Der Wert sollte positiv sein, um die Achsenbeschriftungen sichtbar zu machen.

# Show the plot
plt.suptitle('Bancroft weather parameters histograms', fontsize=24)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])  # Angepasst, um den Titel und die Plots richtig darzustellen
plt.show()
plt.savefig('../output/Bancroft weather parameters histograms.png')

<Figure size 2000x1500 with 0 Axes>


def plot_dotplot(data, parameter, unit, ax):
    sns.stripplot(x=data[parameter], ax=ax, color='blue', size=8, jitter=True)
    ax.set_title(f'Dot Plot for {parameter}')
    ax.set_ylabel('Frequency')
    ax.set_xlabel(unit)

# Set up the subplots with shared y-axis
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=[20,15])

# Plot dot plots for all parameters
for i, (parameter, unit) in enumerate(parameters):
    row, col = divmod(i, 3)
    plot_dotplot(df, parameter, unit, axes[row, col])

# Manually set bottom margin to make x-axis visible
plt.subplots_adjust(bottom=-0.1)

# Show the plot
plt.suptitle('Bancroft weather parameters dot plots', fontsize=24)
plt.tight_layout()
plt.show()
plt.savefig('../output/Bancroft weather parameters dot plots.png')

<Figure size 2000x1500 with 0 Axes>


# Set up the subplots with shared x-axis
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20, 15), constrained_layout=True)

# List of relevant weather parameters
relevant_parameters = ['avg_temp', 'min_wet_bulb_temp', 'avg_dewpoint', 'avg_temp_change', 'avg_windspd', 'max_windgust',
                       'avg_winddir', 'max_cumulative_precip', 'max_snow_density_6', 'max_cumulative_snow',
                       'max_cumulative_ice', 'avg_pressure_change']

sns.set(style="darkgrid")

# Loop through each parameter and create a boxplot
for i, parameter in enumerate(relevant_parameters):
    sns.boxplot(x=df[parameter], ax=axes.flatten()[i], color='blue')
    axes.flatten()[i].set_title(f'Box Plot for {parameter}')
    axes.flatten()[i].set_ylabel(parameter)

# Adjust layout
plt.suptitle('Bancroft weather parameters box plots', fontsize=24)
plt.tight_layout()
plt.show()
plt.savefig('../output/Bancroft weather parameters box plots.png')

<Figure size 2000x1500 with 0 Axes>


def plot_wep_features_distogram(df, label_col):

    f, axes = plt.subplots(4, 3, constrained_layout = True, figsize=[20,15])
    sns.set(style="darkgrid")
    for i in df[label_col].unique():
        df_cluster_distogram = df.loc[df[label_col] == i]
        axes = axes.ravel()
        num_label = df_cluster_distogram[label_col].value_counts()[i]
        for j, feature in enumerate(relevant_parameters):
            sns.distplot(df_cluster_distogram[feature], label=f'{i}: {num_label}', ax=axes[j])
            
    f.suptitle("Bancroft Distribution of Weather Features by Weather Event Profiles in Distograms", fontsize=24)
    
    axes[0].set(title="Avg temperature in Kelvin", xlim=(230, 310), xlabel="")
    h, l = axes[0].get_legend_handles_labels()
    axes[1].set(title="Min wet bulb temperature in Kelvin", xlim=(230, 310), xlabel="")
    axes[2].set(title="Avg dewpoint in Kelvin", xlim=(230, 310), xlabel="")
    axes[3].set(title="Avg temperature change in Kelvin", xlim=(-1, 1), xlabel="")
    axes[4].set(title="Avg wind speed in m/s", xlim=(0, 10), xlabel="")
    axes[5].set(title="Max windgust in m/s", xlim=(0, 30), xlabel="")
    axes[6].set(title="Avg winddirection in degrees", xlim=(0, 360), xlabel="")
    axes[7].set(title="Max cumulative precipitation in mm", xlim=(0, 100), ylim=(0, 0.1), xlabel="")
    axes[8].set(title="Max snow density (when snow exceeds 6.35mm in kg/m^3)", xlim=(0, 800), ylim=(0, 0.01), xlabel="")
    axes[9].set(title="Max cumulative snow accretion in mm", xlim=(0, 1000), ylim=(0, 0.009), xlabel="")
    axes[10].set(title="Max cumulative ice accretion in mm", xlim=(0, 10), ylim=(0, 0.8), xlabel="")
    axes[11].set(title="Avg pressure change in Pa",  xlim=(-200, 200), xlabel="")  
    f.legend(h, l, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=5)
    

data = df.dropna(subset=['wep'], how='any').sort_values('wep')
plot_wep_features_distogram(data, 'wep')
plt.savefig('../output/Bancroft_wep_parameter_distograms.png')


def plot_wep_features_dotplot(df, label_col, relevant_parameters):
    unique_labels = df[label_col].unique()
    palette = sns.color_palette() 
    label_colors = dict(zip(unique_labels, palette))

    f, axes = plt.subplots(4, 3, constrained_layout=True, figsize=[20, 15])
    sns.set(style="darkgrid")

    axes = axes.ravel()
    for i in unique_labels:
        df_cluster_dotplot = df.loc[df[label_col] == i]
        for j, feature in enumerate(relevant_parameters):
            sns.stripplot(x=df_cluster_dotplot[feature], color=label_colors[i], ax=axes[j])

    f.suptitle("Bancroft Distribution of Weather Features by Weather Event Profiles in Dotplots", fontsize=24)
    patches = [plt.Line2D([0], [0], color=label_colors[label], marker='o', linestyle='', markersize=10) for label in unique_labels]
    f.legend(patches, unique_labels, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=len(unique_labels))

    axes[0].set(title="Avg temperature in Kelvin", xlim=(230, 310), xlabel="")
    axes[0].legend(loc="best")
    axes[1].set(title="Min wet bulb temperature in Kelvin", xlim=(230, 310), xlabel="")
    axes[2].set(title="Avg dewpoint in Kelvin", xlim=(230, 310), xlabel="")
    axes[3].set(title="Avg temperature change in Kelvin", xlim=(-1, 1), xlabel="")
    axes[4].set(title="Avg wind speed in m/s", xlim=(0, 10), xlabel="")
    axes[5].set(title="Max windgust in m/s", xlim=(0, 30), xlabel="")
    axes[6].set(title="Avg winddirection in degrees", xlim=(0, 360), xlabel="")
    axes[7].set(title="Max cumulative precipitation in mm", xlim=(0, 100), ylim=(0, 0.1), xlabel="")
    axes[8].set(title="Max snow density (when snow exceeds 6.35mm in kg/m^3)", xlim=(0, 800), ylim=(0, 0.01), xlabel="")
    axes[9].set(title="Max cumulative snow accretion in mm", xlim=(0, 1000), ylim=(0, 0.009), xlabel="")
    axes[10].set(title="Max cumulative ice accretion in mm", xlim=(0, 10), ylim=(0, 0.8), xlabel="")
    axes[11].set(title="Avg pressure change in Pa",  xlim=(-200, 200), xlabel="")  

data = df.dropna(subset=['wep'], how='any').sort_values('wep')
plot_wep_features_dotplot(data, 'wep', relevant_parameters)
plt.savefig('../output/Bancroft_wep_parameter_dotplots.png')

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.


def plot_wep_features_boxplot(df, label_col, relevant_parameters):
    unique_labels = df[label_col].unique()
    palette = sns.color_palette()
    label_colors = dict(zip(unique_labels, palette))

    f, axes = plt.subplots(4, 3, constrained_layout=True, figsize=[20, 15])
    sns.set(style="darkgrid")

    for j, feature in enumerate(relevant_parameters):
        ax = axes.flatten()[j]
        sns.boxplot(x=df[label_col], y=df[feature], ax=ax, palette=label_colors)

        ax.set_title(f'Box Plot for {feature}')
        ax.set_ylabel(feature)
        ax.set_xlabel('')
        ax.set_xticks([])
        ax.set_xticklabels([])
    f.suptitle("Bancroft Distribution of Weather Features by Weather Event Profiles in Boxplots", fontsize=24)
    handles = [plt.Line2D([0], [0], color=label_colors[label], marker='o', linestyle='', markersize=10) for label in unique_labels]
    f.legend(handles, unique_labels, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=len(unique_labels))
    plt.tight_layout()
    plt.show()

data = df.dropna(subset=['wep'], how='any').sort_values('wep')
plot_wep_features_boxplot(data, 'wep', relevant_parameters)
plt.savefig('../output/Bancroft_wep_parameter_boxplots.png')

<Figure size 2000x1500 with 0 Axes>


# If 'Blue sky day' is present in 'wep', use 'label0'; otherwise, use 'wep'
df['display_label'] = df.apply(lambda row: f"{row['wep']}: '{row['label0']}'" if row['wep'] == 'Blue sky day' else row['wep'], axis=1)

# Count occurrences of weather events
weather_counts = df['display_label'].value_counts()

# Count occurrences of weather events excluding 'Blue sky day'
not_blue_sky_day_counts = df[df['wep'] != 'Blue sky day']['wep'].value_counts()

# Determine events to be exploded
explode_events = set(not_blue_sky_day_counts.index)

# Create explode list for the first chart
explode = [0.1 if event in explode_events else 0 for event in weather_counts.index]

# Total number of all samples
total_samples = len(df)

# Create the plot
plt.figure(figsize=(16, 8))

# Group categories and adjust labels
lower_half_threshold = weather_counts.quantile(0.5)
grouped_labels = [label if weather_counts[label] > lower_half_threshold else 'Other' for label in weather_counts.index]

# Pie chart for all weather events
plt.subplot(1, 2, 1)
patches, texts, autotexts = plt.pie(weather_counts, labels=grouped_labels, autopct=lambda p: f'{p:.1f}%\n({int(p/100*total_samples)})', startangle=90, explode=explode)
plt.title(f'Distribution of All Weather Events\nTotal Number of Observations: {total_samples}')

# Adjust text size and hide overlapping labels
for text, autotext in zip(texts, autotexts):
    if text.get_text() == 'Other':
        text.set_visible(False)
    text.set_fontsize(8)
    autotext.set_fontsize(8)

# Enlarged view for 'wep' != 'Blue sky day'
plt.subplot(1, 2, 2)
total_not_blue_samples = not_blue_sky_day_counts.sum()
plt.pie(not_blue_sky_day_counts, labels=not_blue_sky_day_counts.index, autopct=lambda p: f'{p:.1f}%\n({int(p/100*total_not_blue_samples)})', startangle=90)
plt.title(f'Distribution of Extreme Weather Events\nTotal Number of Non-Blue-Sky Observations: {total_not_blue_samples}')
plt.show()
plt.savefig('../output/Distribution of weather events pie chart.png')

<Figure size 2000x1500 with 0 Axes>


df_datetime = df_full.copy()

# Grouping the data by year and month and calculating the average in Celsius
grouped_monthly_avg_celsius = df_full.groupby([df_full['run_datetime'].dt.year, df_full['run_datetime'].dt.month])['avg_temp_celsius'].mean()

# Conversion of the grouped data into a format suitable for the radar type
radar_data_celsius = grouped_monthly_avg_celsius.unstack(level=1)

# Selecting the color chart
cmap = plt.cm.coolwarm

# Determining the temperatures in July and scaling for the color chart
july_temps = radar_data_celsius.loc[:, 7]  # 7 steht für Juli
min_july_temp = july_temps.min()
max_july_temp = july_temps.max()
scaled_july_temps = (july_temps - min_july_temp) / (max_july_temp - min_july_temp)
colors = [cmap(temp) for temp in scaled_july_temps]  # Liste der Farben

# Creating the radar chart with Celsius data
angles_celsius = [n / 12 * 2 * math.pi for n in range(12)]  # 12 Monate
angles_celsius += angles_celsius[:1]  # Schließt den Kreis

# Month names
month_names = ["Jan", "Feb", "Mär", "Apr", "Mai", "Jun", "Jul", "Aug", "Sep", "Okt", "Nov", "Dez"]

# Creating the radar chart
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(polar=True))
for i, year in enumerate(radar_data_celsius.index):
    values_celsius = radar_data_celsius.loc[year].tolist()
    values_celsius += values_celsius[:1]  # Closing the circle
    ax.plot(angles_celsius, values_celsius, linewidth=1, linestyle='solid', label=f'Year {year}', color=colors[i])

# Add month names as labels
plt.xticks(angles_celsius[:-1], month_names, color='grey', size=12)

plt.title('Temperature change over the months of all years in the dataset (in Celsius)', size=20, y=1.1)
plt.legend(loc='upper right', bbox_to_anchor=(0.1, 0.1))
plt.show()
plt.savefig('../output/Temperature change radar chart.png')

<Figure size 2000x1500 with 0 Axes>


weather_counts = df['wep'].value_counts()
print(weather_counts)

wep
Blue sky day                                           52633
Mild snowfall                                           4497
Moderate snowfall                                       2920
Moderate rain                                           2630
Heavy snowfall with accumulated snow                    1040
Frondurchlauf / Continuous freezing rain                 934
Storm with freezing rain / Heavy snow- and icestorm      430
Snowstorm with high precipitation                        261
Name: count, dtype: int64


# Convert the 'run_datetime' column to a datetime object
df['run_datetime'] = pd.to_datetime(df['run_datetime'])

# Aggregate the average wind speed and temperature for each month across all years
avg_by_month = df.groupby(df['run_datetime'].dt.to_period('M')).agg({
    'avg_windspd': 'mean',
    'avg_temp': 'mean'
}).reset_index()

# Ensure the presence of a date column
avg_by_month['Month'] = avg_by_month['run_datetime'].dt.strftime('%b')

# Adjust the order of the months
months_order = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
avg_by_month['Month'] = pd.Categorical(avg_by_month['Month'], categories=months_order, ordered=True)
avg_by_month = avg_by_month.sort_values('Month')

# Plot as a Line Chart
fig, ax1 = plt.subplots(figsize=(15, 8))

# Line for wind speed on the left y-axis
color = 'tab:blue'
ax1.set_xlabel('Month')
ax1.set_ylabel('Average windspeed (m/s)', color=color)
ax1 = sns.lineplot(x='Month', y='avg_windspd', data=avg_by_month, sort=False, label='Wind Speed', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# Second y-axis for temperature
ax2 = ax1.twinx()
color = 'tab:red'
ax2.set_ylabel('Average temperature (°C)', color=color)
ax2 = sns.lineplot(x='Month', y='avg_temp', data=avg_by_month, sort=False, label='Temperature', color=color)
ax2.tick_params(axis='y', labelcolor=color)

plt.title('Average Wind Speed and Temperature per Month Across All Years')

# Adding the legend manually
lines_1, labels_1 = ax1.get_legend_handles_labels()
lines_2, labels_2 = ax2.get_legend_handles_labels()
plt.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper right')

plt.show()
plt.savefig('../output/verage wind speed and temperature per month across all years.png')

<Figure size 2000x1500 with 0 Axes>


# Display a scatter plot in a diagram with subplots for all parameter pairs
fig, axes = plt.subplots(len(relevant_parameters), len(relevant_parameters), figsize=(50, 50), gridspec_kw={'hspace': 0.8, 'wspace': 0.4})

# Function for calculating the correlation coefficient and displaying the scatterplot
def plot_scatter_and_correlation(ax, x, y, data):
    sns.scatterplot(x=x, y=y, data=data, ax=ax)
    correlation_coefficient = data[x].corr(data[y])
    ax.set_title(f'Scatterplot of {x} and {y}', fontsize=16)
    ax.set_xlabel(x, fontsize=14)
    ax.set_ylabel(y, fontsize=14)
    ax.text(0.5, -0.27, f'Correlation Coefficient: {correlation_coefficient:.2f}', transform=ax.transAxes, ha='center', fontsize=12)

# Loop over all parameters 
for i, feature_x in enumerate(relevant_parameters):
    for j, feature_y in enumerate(relevant_parameters):
        if i != j:  # Damit die Diagonale übersprungen wird (x gegen x)
            plot_scatter_and_correlation(axes[i, j], feature_x, feature_y, df)

plt.tight_layout()
plt.show()
plt.savefig('../output/Bancroft Association Plots.png')

<Figure size 2000x1500 with 0 Axes>


def plot_wind_speed_and_avg_temp_by_direction(df_full, df, wind_direction_label, avg_temp_label, avg_windspd_label):
    # Calculate the average temperature for each wind direction
    avg_temp_by_direction = df_full.groupby(wind_direction_label)[avg_temp_label].mean().reset_index()

    # Order cardinal directions clockwise
    directions_ordered = ['North', 'Northeast', 'East', 'Southeast', 'South', 'Southwest', 'West', 'Northwest']
    avg_temp_by_direction[wind_direction_label] = pd.Categorical(avg_temp_by_direction[wind_direction_label], categories=directions_ordered, ordered=True)
    avg_temp_by_direction = avg_temp_by_direction.sort_values(wind_direction_label)

    # Color scale for temperature
    color_scale = plt.cm.get_cmap('coolwarm')

    # Create the bar chart
    fig, ax = plt.subplots(figsize=(8, 6))

    # Loop through the data and create bars with corresponding color
    for i, row in avg_temp_by_direction.iterrows():
        color = color_scale(row[avg_temp_label] / avg_temp_by_direction[avg_temp_label].max())
        wind_direction = row[wind_direction_label]
        average_wind_speed = df[df[wind_direction_label] == wind_direction][avg_windspd_label].mean()
        ax.bar(wind_direction, average_wind_speed, color=color, edgecolor='black')

    # Chart settings
    ax.set_xlabel('Wind Direction')
    ax.set_ylabel('Wind Speed (m/s)')
    ax.set_title('Wind Speeds and Average Temperatures by Wind Direction')

    # Create a separate Axes for the colorbar
    cax = fig.add_axes([0.95, 0.1, 0.04, 0.79])

    # Add color legend
    temperature_min = avg_temp_by_direction[avg_temp_label].min()
    temperature_max = avg_temp_by_direction[avg_temp_label].max()
    sm = plt.cm.ScalarMappable(cmap=color_scale, norm=plt.Normalize(vmin=temperature_min, vmax=temperature_max))
    sm.set_array([])  # empty data for the color legend
    cbar = plt.colorbar(sm, cax=cax, orientation='vertical')
    cbar.set_label('Average Temperature (°C)')

    # Show the chart
    plt.show()

# Beispielaufruf der Funktion
plot_wind_speed_and_avg_temp_by_direction(df_full, df, 'wind_direction_label', 'avg_temp_celsius', 'avg_windspd')
plt.savefig('../output/Bancroft Wind Speeds and Average Temperatures by Wind Direction.png')

<Figure size 2000x1500 with 0 Axes>


import plotly.express as px
from sklearn.decomposition import PCA

feature_columns = df.columns
columns_to_exclude = ['wep', 'run_datetime', 'avg_temp_celsius', 'avg_winddir', 'wind_direction_label', 'label0', 'label1', 'label2', 'year', 'month', 'display_label']
feature_columns = [col for col in feature_columns if col not in columns_to_exclude]

X = df[feature_columns].dropna()
X_standardized = StandardScaler().fit_transform(X)

# Apply PCA to reduce dimensionality to 3 components
pca = PCA(n_components=3)
X_pca = pca.fit_transform(X_standardized)

# Create a DataFrame with PCA results and 'wep' label
pca_df = pd.DataFrame(data=X_pca, columns=['PC1', 'PC2', 'PC3'])
pca_df['wep'] = df.loc[X.index, 'wep']

# Create an interactive 3D scatter plot using Plotly Express with smaller dots
fig = px.scatter_3d(pca_df, x='PC1', y='PC2', z='PC3', color='wep', symbol='wep', size_max=2, opacity=0.7,
                    labels={'wep': 'Weather'}, title='Interactive 3D PCA Plot of Weather Data')

# Update layout for better visualization
fig.update_layout(margin=dict(l=0, r=0, b=0, t=0))

# Update marker size for better visibility
fig.update_traces(marker=dict(size=2),
                  selector=dict(mode='markers'))

# Increase marker size in the legend
fig.update_layout(legend=dict(
    itemsizing='constant'
))

# Show the plot
fig.show()


#import dependencies for models (import2)
from sklearn import svm

from datetime import datetime

from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error, classification_report, confusion_matrix, accuracy_score, precision_score, recall_score, f1_score

from sklearn.linear_model import LinearRegression, LogisticRegression, SGDRegressor

from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier, ExtraTreesClassifier, RandomForestClassifier, GradientBoostingClassifier

from sklearn.tree import DecisionTreeClassifier

from sklearn.neighbors import KNeighborsClassifier

from sklearn.svm import SVR

from sklearn.pipeline import make_pipeline

from xgboost import XGBClassifier
from lightgbm import LGBMClassifier

import statsmodels.api as sm

import pickle


# Create Linear Regression model
lr_model_1 = LinearRegression()

# Create Gradient Boosting Regressor
gb_model_1 = GradientBoostingRegressor(n_estimators=100, random_state=42)

# Create Stochastic Gradient Descent Regressor
sgd_model_1 = SGDRegressor(max_iter=1000, random_state=42)

# Create Support Vector Regressor
svr_model_1 = SVR()


# Select the features (independent variable)
X = df['avg_windspd'].values.reshape(-1, 1)

# Select the dependent variable
y = df['avg_temp']

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Print the shapes of the sets
print("X_train shape:", X_train.shape)
print("X_test shape:", X_test.shape)
print("y_train shape:", y_train.shape)
print("y_test shape:", y_test.shape)

X_train shape: (52276, 1)
X_test shape: (13069, 1)
y_train shape: (52276,)
y_test shape: (13069,)


# fit Linear Regression model
lr_model_fit_1 = lr_model_1.fit(X_train, y_train)

# fit Gradient Boosting Regressor
gb_model_fit_1 = gb_model_1.fit(X_train, y_train)

# fit Stochastic Gradient Descent Regressor
sgd_model_fit_1 = sgd_model_1.fit(X_train, y_train)

# fit Support Vector Regressor
svr_model_fit_1 = svr_model_1.fit(X_train, y_train)


# Function to evaluate and print model performance
def evaluate_model(model_name, model, X_test, y_test):
    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    mae = mean_absolute_error(y_test, y_pred)
    print(f'{model_name} Model:')
    print(f'Mean Squared Error: {mse:.2f}')
    print(f'Mean Absolute Error: {mae:.2f}')
    print()

# Evaluate Linear Regression model
evaluate_model('Linear Regression', lr_model_fit_1, X_test, y_test)

# Evaluate Gradient Boosting model
evaluate_model('Gradient Boosting', gb_model_fit_1, X_test, y_test)

# Evaluate Stochastic Gradient Descent model
evaluate_model('Stochastic Gradient Descent', sgd_model_fit_1, X_test, y_test)

# Evaluate Support Vector Regression model
evaluate_model('Support Vector Regression', svr_model_fit_1, X_test, y_test)

Linear Regression Model:
Mean Squared Error: 127.78
Mean Absolute Error: 9.63

Gradient Boosting Model:
Mean Squared Error: 127.72
Mean Absolute Error: 9.63

Stochastic Gradient Descent Model:
Mean Squared Error: 127.78
Mean Absolute Error: 9.63

Support Vector Regression Model:
Mean Squared Error: 129.85
Mean Absolute Error: 9.57


from scipy.stats import zscore

# Function to create subplots for model visualization and outlier analysis
def visualize_model(model_name, model, X_test, y_test, subplot_idx):
    y_pred = model.predict(X_test)

    # Calculate residuals
    residuals = y_test - y_pred

    # Calculate MSE and MAE
    mse = mean_squared_error(y_test, y_pred)
    mae = mean_absolute_error(y_test, y_pred)

    # Detect outliers using Z-scores
    z_scores = zscore(residuals)
    outliers = np.abs(z_scores) > 2  # Threshold for considering a point as an outlier

    # Create subplots for actual vs predicted values
    plt.subplot(2, 4, subplot_idx)
    plt.scatter(X_test[outliers], y_test[outliers], color='orange', label='Outliers')
    plt.scatter(X_test[~outliers], y_test[~outliers], color='blue', label='Actual')
    plt.scatter(X_test, y_pred, color='red', label='Predicted', alpha=0.7)
    plt.plot(X_test, y_pred, color='green', linewidth=2)
    plt.title(f'{model_name} Model\nMSE: {mse:.2f}, MAE: {mae:.2f}')
    plt.xlabel('Average Wind Speed')
    plt.ylabel('Average Temperature')

    # Create subplots for residual plots
    plt.subplot(2, 4, subplot_idx + 4)
    plt.scatter(X_test[outliers], residuals[outliers], color='orange', label='Outliers')
    plt.scatter(X_test[~outliers], residuals[~outliers], color='purple')
    plt.axhline(y=0, color='black', linestyle='--', linewidth=2)
    plt.title(f'Residuals\nMean Residual: {np.mean(residuals):.2f}')
    plt.xlabel('Average Wind Speed')
    plt.ylabel('Residuals')

# Create subplots
plt.figure(figsize=(16, 8))

# Visualize Linear Regression model with outliers
visualize_model('Linear Regression', lr_model_fit_1, X_test, y_test, 1)

# Visualize Gradient Boosting model with outliers
visualize_model('Gradient Boosting', gb_model_fit_1, X_test, y_test, 2)

# Visualize Stochastic Gradient Descent model with outliers
visualize_model('Stochastic Gradient Descent', sgd_model_fit_1, X_test, y_test, 3)

# Visualize Support Vector Regression model with outliers
visualize_model('Support Vector Regression', svr_model_fit_1, X_test, y_test, 4)

# Adjust layout for better visualization
plt.tight_layout()

# Create a single legend at the top of the entire plot
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.1), fancybox=True, shadow=True, ncol=4)

plt.show()
plt.savefig('../output/Linear Regression Analysis Temperature and Wind Modeling Results.png')

<Figure size 2000x1500 with 0 Axes>


# Create a meaningful name with a timestamp
model_name_lr_1 = f"lr_model_fit_1_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_gb_1 = f"gb_model_fit_1_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_sgd_1 = f"sgd_model_fit_1_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_svr_1 = f"svr_model_fit_1_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name_lr_1, 'wb') as f:
    pickle.dump(lr_model_fit_1, f)
with open(folder_path + model_name_gb_1, 'wb') as f:
    pickle.dump(gb_model_fit_1, f)
with open(folder_path + model_name_sgd_1, 'wb') as f:
    pickle.dump(sgd_model_fit_1, f)
with open(folder_path + model_name_svr_1, 'wb') as f:
    pickle.dump(svr_model_fit_1, f)


#drop columns with datatype object
df_feature_selection = df.drop(columns=['run_datetime','wep','wind_direction_label','label2','display_label'])

# Calculate the correlation matrix
correlation_matrix = df_feature_selection.corr()

# Create a heatmap for the correlation matrix
plt.figure(figsize=(20, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('Correlation Matrix of Variables')
plt.show()
plt.savefig('../output/Linear Regression Analysis Multiple Predictors Correlation Matrix of Variables.png')

<Figure size 2000x1500 with 0 Axes>


#set "avg_temp_change" to 0 in the first row, as there is no temp change
df_feature_selection.at[0, 'avg_temp_change'] = 0

#drop columns, which are not numeric and not needed
df_feature_selection = df_feature_selection.drop(columns = ['min_wet_bulb_temp', 'avg_temp_change', 'max_windgust' , 'avg_dewpoint' , 'avg_temp_celsius','avg_pressure_change','label0'] )

#show adjusted df
df_feature_selection.head()


# Extracting the target variable and predictor variables
y = df_feature_selection['label1']
#drop target variable and strongly correlated varaibles
X = df_feature_selection.drop(columns=['label1'])

# Initializing an empty list to store the selected variables
selected_vars = []
# Initializing an empty list to store adjusted R-squared scores for each step
adjusted_r2_scores = []

# Performing forward selection
for i in range(len(X.columns)):
    best_score = -1  # Initialize with a value that will be replaced on the first iteration
    best_var = None

    # Iterate through remaining variables
    for var in X.columns:
        if var not in selected_vars:
            # Include the variable in the model
            current_vars = selected_vars + [var]
            X_current = X[current_vars]
            X_current = sm.add_constant(X_current)  # Adding a constant for the intercept

            # Fit the model
            model = sm.OLS(y, X_current).fit()
            current_score = model.rsquared_adj

            # Update the best variable if the score is higher
            if current_score > best_score:
                best_score = current_score
                best_var = var

    # Update selected variables and store results
    selected_vars.append(best_var)
    adjusted_r2_scores.append(best_score)
    print(f"Step {i + 1}: Added {best_var}, Adjusted R-squared: {best_score:.4f}")

# Display the results
results_df = pd.DataFrame({'Variable': selected_vars, 'Adjusted R-squared': adjusted_r2_scores})
print("\nFinal Forward Selection Results:")
print(results_df)

Step 1: Added max_snow_density_6, Adjusted R-squared: 0.4062
Step 2: Added avg_temp, Adjusted R-squared: 0.5109
Step 3: Added max_cumulative_snow, Adjusted R-squared: 0.5506
Step 4: Added max_cumulative_ice, Adjusted R-squared: 0.5527
Step 5: Added avg_windspd, Adjusted R-squared: 0.5541
Step 6: Added avg_winddir, Adjusted R-squared: 0.5551
Step 7: Added max_cumulative_precip, Adjusted R-squared: 0.5556
Step 8: Added avg_winddir_sin, Adjusted R-squared: 0.5556
Step 9: Added avg_winddir_cos, Adjusted R-squared: 0.5556

Final Forward Selection Results:
                Variable  Adjusted R-squared
0     max_snow_density_6            0.406223
1               avg_temp            0.510857
2    max_cumulative_snow            0.550552
3     max_cumulative_ice            0.552674
4            avg_windspd            0.554109
5            avg_winddir            0.555149
6  max_cumulative_precip            0.555629
7        avg_winddir_sin            0.555632
8        avg_winddir_cos            0.555626


X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Train a model with all features
model = RandomForestClassifier(random_state=42)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
initial_accuracy = accuracy_score(y_test, y_pred)

# Perform Backward Feature Elimination
all_features = list(X.columns)
selected_features = all_features.copy()

while len(selected_features) > 1:  # Stopping criterion: keep at least one feature
    # Train model with the current set of features
    model.fit(X_train[selected_features], y_train)
    y_pred = model.predict(X_test[selected_features])
    
    # Evaluate model performance
    current_accuracy = accuracy_score(y_test, y_pred)
    
    # If removing a feature improves performance, remove the least important one
    if current_accuracy >= initial_accuracy:
        break
    else:
        least_important_feature = min(selected_features, key=lambda f: model.feature_importances_[selected_features.index(f)])
        selected_features.remove(least_important_feature)

print("Optimal Feature-Subset:")
for feature in selected_features:
    print(feature)
print("Accuracy: ",current_accuracy)

Optimal Feature-Subset:
avg_temp
avg_windspd
avg_winddir
avg_winddir_sin
avg_winddir_cos
max_cumulative_precip
max_snow_density_6
max_cumulative_snow
max_cumulative_ice
Accuracy:  0.9996939322059836


# Extract relevant columns from the DataFrame
results_df = df[["run_datetime","avg_temp", "avg_winddir", "avg_windspd","avg_winddir_sin","avg_winddir_cos","max_cumulative_precip","max_snow_density_6","max_cumulative_snow","max_cumulative_ice"]]

# Set the 'run_datetime' column as the DataFrame index and convert it to datetime format
results_df['run_datetime'] = pd.to_datetime(results_df['run_datetime'])
results_df.set_index('run_datetime', inplace=True)

# Remove rows with duplicated indices, keeping only the first occurrence
results_df = results_df.loc[~results_df.index.duplicated(keep='first')]

# Set "avg_temp_change" to 0 in the first row, as there is no temp change
results_df.iloc[0, results_df.columns.get_loc("avg_temp")] = 0

# Display the first few rows of the resulting DataFrame
results_df.head()


# Extracting the target variable and predictor variables
y = results_df['avg_temp']
X = results_df[['avg_windspd','avg_winddir']]

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


#Training (Multiple Linear Regression)
regr = LinearRegression()
regr.fit(X, y)

# Fit Model (Multiple Linear Regression)
y_pred = regr.predict(X_test)

#Evaluate on Test Set (Multiple Linear Regression)
mse = mean_squared_error(y_test, y_pred)
print('Mean Squared Error:', mse)

mae = mean_absolute_error(y_test, y_pred)
print('Mean Absolute Error:', mae)

r2 = r2_score(y_test, y_pred)
print('R-squared:', r2)

Mean Squared Error: 128.267893719392
Mean Absolute Error: 9.658907432391409
R-squared: 0.015238184007991595


# Create subplots
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Scatter Plot
axes[0].scatter(y_test, y_pred)
axes[0].set_xlabel('True Values')
axes[0].set_ylabel('Predictions')
axes[0].set_title('Scatter Plot of True vs. Predicted Values')

# Residual Plot
residuals = y_test - y_pred
axes[1].scatter(y_pred, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_xlabel('Predicted Values')
axes[1].set_ylabel('Residuals')
axes[1].set_title('Residual Plot')

# Histogram of Residuals
sns.histplot(residuals, kde=True, ax=axes[2])
axes[2].set_xlabel('Residuals')
axes[2].set_ylabel('Frequency')
axes[2].set_title('Histogram of Residuals')

# Add title
fig.suptitle('MLR with target variable "Average Temperature" and predictor variables "Average Windspeed" and "Wind Direction"')

# Adjust layout
plt.tight_layout()

# Show plots
plt.show()


# Create a meaningful name with a timestamp
model_name = f"MLR_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name, 'wb') as f:
    pickle.dump(regr, f)

print(f"The MLR model has been saved to '{folder_path + model_name}'.")

The MLR model has been saved to '../models/MLR_model_20231229_1802.pkl'.


# Extracting the target variable and predictor variables
y = results_df['avg_temp']
X = results_df.drop(columns=['avg_temp'])

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


#Training (Multiple Linear Regression)
regr = LinearRegression()
regr.fit(X, y)

# Fit Model (Multiple Linear Regression)
y_pred = regr.predict(X_test)

#Evaluate on Test Set (Multiple Linear Regression)
mse = mean_squared_error(y_test, y_pred)
print('Mean Squared Error:', mse)

mae = mean_absolute_error(y_test, y_pred)
print('Mean Absolute Error:', mae)

r2 = r2_score(y_test, y_pred)
print('R-squared:', r2)

Mean Squared Error: 103.25403414868913
Mean Absolute Error: 8.390307032430659
R-squared: 0.20727917775583216


# Create subplots
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Scatter Plot
axes[0].scatter(y_test, y_pred)
axes[0].set_xlabel('True Values')
axes[0].set_ylabel('Predictions')
axes[0].set_title('Scatter Plot of True vs. Predicted Values')

# Residual Plot
residuals = y_test - y_pred
axes[1].scatter(y_pred, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_xlabel('Predicted Values')
axes[1].set_ylabel('Residuals')
axes[1].set_title('Residual Plot')

# Histogram of Residuals
sns.histplot(residuals, kde=True, ax=axes[2])
axes[2].set_xlabel('Residuals')
axes[2].set_ylabel('Frequency')
axes[2].set_title('Histogram of Residuals')

# Add title
fig.suptitle('MLR with target variable "Average Temperature" and multiple predictor variables')

# Adjust layout
plt.tight_layout()

# Show plots
plt.show()


# Create a meaningful name with a timestamp
model_name = f"MLR_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name, 'wb') as f:
    pickle.dump(regr, f)

print(f"The MLR model has been saved to '{folder_path + model_name}'.")

The MLR model has been saved to '../models/MLR_model_20231229_1802.pkl'.


# Extracting the target variable and predictor variables
y = results_df['avg_windspd']
X = results_df.drop(columns=['avg_windspd'])

# Splitting the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


#Training (Multiple Linear Regression)
regr = LinearRegression()
regr.fit(X, y)

#  Fit Model (Multiple Linear Regression)
y_pred = regr.predict(X_test)

#Evaluate on Test Set (Multiple Linear Regression)
mse = mean_squared_error(y_test, y_pred)
print('Mean Squared Error:', mse)

mae = mean_absolute_error(y_test, y_pred)
print('Mean Absolute Error:', mae)

r2 = r2_score(y_test, y_pred)
print('R-squared:', r2)

Mean Squared Error: 0.5038677092424138
Mean Absolute Error: 0.5636674787101482
R-squared: 0.1563048200931636


# Create subplots
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Scatter Plot
axes[0].scatter(y_test, y_pred)
axes[0].set_xlabel('True Values')
axes[0].set_ylabel('Predictions')
axes[0].set_title('Scatter Plot of True vs. Predicted Values')

# Residual Plot
residuals = y_test - y_pred
axes[1].scatter(y_pred, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_xlabel('Predicted Values')
axes[1].set_ylabel('Residuals')
axes[1].set_title('Residual Plot')

# Histogram of Residuals
sns.histplot(residuals, kde=True, ax=axes[2])
axes[2].set_xlabel('Residuals')
axes[2].set_ylabel('Frequency')
axes[2].set_title('Histogram of Residuals')

# Add title
fig.suptitle('Multiple Linear Regression (MLR) with target variable "Average Winddirection" and multiple predictor variables')

# Adjust layout
plt.tight_layout()

# Show plots
plt.show()


# Create a meaningful name with a timestamp
model_name = f"MLR_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name, 'wb') as f:
    pickle.dump(regr, f)

print(f"The MLR model has been saved to '{folder_path + model_name}'.")

The MLR model has been saved to '../models/MLR_model_20231229_1802.pkl'.


# Extract relevant columns from the DataFrame
df_temp_wind = df_full[["run_datetime", "avg_temp", "avg_windspd", "avg_winddir", "avg_windgust"]]

# Set the 'run_datetime' column as the DataFrame index and convert it to datetime format
df_temp_wind.index = pd.to_datetime(df_temp_wind['run_datetime'])

# Drop the original 'run_datetime' column as it's now set as the index
df_temp_wind.drop(columns='run_datetime', inplace=True)

# Remove rows with duplicated indices, keeping only the first occurrence
df_temp_wind = df_temp_wind.loc[~df_temp_wind.index.duplicated(keep='first')]

# Display the first few rows of the resulting DataFrame
df_temp_wind.head()


from statsmodels.tsa.seasonal import seasonal_decompose

def decompose_and_plot(data, variable_name, position):
    # Decompose the time series data
    decompose_data = seasonal_decompose(data.asfreq('M'), model="additive", extrapolate_trend='freq')
    
    # Plot the original data
    plt.subplot(4, 4, position)
    data.plot()
    plt.title(f'Original Data of {variable_name}')
    
    # Plot the trend component
    plt.subplot(4, 4, position + 1)
    decompose_data.trend.plot()
    plt.title(f'Trend of {variable_name}')
    
    # Plot the seasonal component
    plt.subplot(4, 4, position + 2)
    decompose_data.seasonal.plot()
    plt.title(f'Seasonal of {variable_name}')
    
    # Plot the residual component
    plt.subplot(4, 4, position + 3)
    decompose_data.resid.plot()
    plt.title(f'Residual of {variable_name}')


# Create a figure with a size of 16x12 inches
plt.figure(figsize=(16, 12))

# Decompose and plot for avg_temp
decompose_and_plot(df_temp_wind["avg_temp"], "avg_temp", 1)

# Decompose and plot for avg_winddir
decompose_and_plot(df_temp_wind["avg_winddir"], "avg_winddir", 5)

# Decompose and plot for avg_windspd
decompose_and_plot(df_temp_wind["avg_windspd"], "avg_windspd", 9)

# Decompose and plot for avg_windgust
decompose_and_plot(df_temp_wind["avg_windgust"], "avg_windgust", 13)

# Adjust layout to prevent overlap
plt.tight_layout()

# Display the plots
plt.show()
plt.savefig('../output/Linear Regression Analysis Multiple Predictors Seasonality and Trend.png')

<Figure size 2000x1500 with 0 Axes>


from statsmodels.tsa.stattools import adfuller

def adf_test(variable, name):
    # Perform the Augmented Dickey-Fuller test
    dftest = adfuller(variable, autolag='AIC')
    
    # Print the results
    print(f"{name}: ")
    print(f"1. ADF : {dftest[0]}")
    print(f"2. P-Value : {dftest[1]}")
    print(f"3. Num Of Lags : {dftest[2]}")
    print(f"4. Num Of Observations Used For ADF Regression and Critical Values Calculation : {dftest[3]}")
    print("5. Critical Values :")
    for key, val in dftest[4].items():
        print(f"\t{key} : {val}")
    print("\n")

# Perform ADF test for each variable
adf_test(df_temp_wind["avg_temp"], "Average Temperature")
adf_test(df_temp_wind["avg_winddir"], "Average Wind Direction")
adf_test(df_temp_wind["avg_windspd"], "Average Wind Speed")
adf_test(df_temp_wind["avg_windgust"], "Average Wind Gust")

Average Temperature: 
1. ADF : -5.6633939102963415
2. P-Value : 9.268984051353608e-07
3. Num Of Lags : 61
4. Num Of Observations Used For ADF Regression and Critical Values Calculation : 65258
5. Critical Values :
	1% : -3.4304502108131385
	5% : -2.8615842913495717
	10% : -2.566793574780786


Average Wind Direction: 
1. ADF : -26.935153074526742
2. P-Value : 0.0
3. Num Of Lags : 54
4. Num Of Observations Used For ADF Regression and Critical Values Calculation : 65265
5. Critical Values :
	1% : -3.4304502000646013
	5% : -2.861584286598995
	10% : -2.5667935722522013


Average Wind Speed: 
1. ADF : -26.388919718643297
2. P-Value : 0.0
3. Num Of Lags : 53
4. Num Of Observations Used For ADF Regression and Critical Values Calculation : 65266
5. Critical Values :
	1% : -3.4304501985292846
	5% : -2.861584285920424
	10% : -2.566793571891019


Average Wind Gust: 
1. ADF : -24.04907190849269
2. P-Value : 0.0
3. Num Of Lags : 61
4. Num Of Observations Used For ADF Regression and Critical Values Calculation : 65258
5. Critical Values :
	1% : -3.4304502108131385
	5% : -2.8615842913495717
	10% : -2.566793574780786


import itertools

# Define response variable
y = df_temp_wind['avg_temp']

# Define predictor variables
predictor_vars = ['avg_windspd', 'avg_winddir', 'avg_windgust']

# Initialize variables to store best AIC and best predictor combination
best_aic = float('inf')
best_predictors = None

# Iterate through all combinations of predictor variables
for r in range(1, len(predictor_vars) + 1):
    for combo in itertools.combinations(predictor_vars, r):
        # Select the current combination of predictors
        current_predictors = list(combo)
        
        # Define the predictors and add a constant
        x = df_temp_wind[current_predictors]
        x = sm.add_constant(x)

        # Fit the regression model
        model = sm.OLS(y, x).fit()

        # Update best predictors if the current AIC is lower
        if model.aic < best_aic:
            best_aic = model.aic
            best_predictors = current_predictors

# Print the best combination of predictors and its AIC
print(f"Best Predictors: {best_predictors}")
print(f"Best AIC: {best_aic}")

Best Predictors: ['avg_windspd', 'avg_winddir']
Best AIC: 502202.826219331


# Define response variable
y = df_temp_wind['avg_temp']

# Define predictor variables
predictor_vars = ['avg_windspd', 'avg_winddir','avg_windgust']

# Initialize variables to store best BIC and best predictor combination
best_bic = float('inf')
best_predictors_bic = None

# Iterate through all combinations of predictor variables
for r in range(1, len(predictor_vars) + 1):
    for combo in itertools.combinations(predictor_vars, r):
        # Select the current combination of predictors
        current_predictors = list(combo)
        
        # Define the predictors and add a constant
        x = df_temp_wind[current_predictors]
        x = sm.add_constant(x)

        # Fit the regression model
        model = sm.OLS(y, x).fit()

        # Update best predictors if the current BIC is lower
        if model.bic < best_bic:
            best_bic = model.bic
            best_predictors_bic = current_predictors

# Print the best combination of predictors and its BIC
print(f"Best Predictors (BIC): {best_predictors_bic}")
print(f"Best BIC: {best_bic}")

Best Predictors (BIC): ['avg_windspd', 'avg_winddir']
Best BIC: 502230.0873799723


from pmdarima import auto_arima

auto_arima_results= auto_arima(df_temp_wind["avg_temp"], 
           start_p=1,
           start_q=1,
           max_p=3,
           max_q=3, # maximum p and q
           m=24,              # frequency of series
           d=None,           # let model determine 'd'
           seasonal=False,   # No Seasonality
           start_P=0, 
           D=1,  
           trace=True,
           error_action='ignore',  
           suppress_warnings=True, 
           stepwise=True)

auto_arima_results.summary()

Performing stepwise search to minimize aic
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=-300996.225, Time=7.82 sec
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=-30853.412, Time=2.62 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=-272301.363, Time=3.90 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=-116812.636, Time=10.22 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=-30855.181, Time=1.35 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=-315512.530, Time=3.67 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=-315504.206, Time=4.01 sec
 ARIMA(3,1,1)(0,0,0)[0] intercept   : AIC=-315506.573, Time=3.23 sec
 ARIMA(2,1,2)(0,0,0)[0] intercept   : AIC=-315513.680, Time=8.53 sec
 ARIMA(1,1,2)(0,0,0)[0] intercept   : AIC=-309847.293, Time=22.26 sec
 ARIMA(3,1,2)(0,0,0)[0] intercept   : AIC=-309557.825, Time=40.41 sec
 ARIMA(2,1,3)(0,0,0)[0] intercept   : AIC=-315534.509, Time=12.44 sec
 ARIMA(1,1,3)(0,0,0)[0] intercept   : AIC=-312901.486, Time=18.11 sec
 ARIMA(3,1,3)(0,0,0)[0] intercept   : AIC=-315657.922, Time=10.30 sec
 ARIMA(3,1,3)(0,0,0)[0]             : AIC=-315828.456, Time=23.75 sec
 ARIMA(2,1,3)(0,0,0)[0]             : AIC=-315536.495, Time=5.94 sec
 ARIMA(3,1,2)(0,0,0)[0]             : AIC=inf, Time=30.87 sec
 ARIMA(2,1,2)(0,0,0)[0]             : AIC=-315515.665, Time=4.10 sec

Best model:  ARIMA(3,1,3)(0,0,0)[0]          
Total fit time: 213.543 seconds


auto_arima_results.plot_diagnostics(figsize=(10,8))
plt.show()


from statsmodels.tsa.statespace.sarimax import SARIMAX

#manually define train and test data
train_df = df_temp_wind[df_temp_wind.index <= datetime(2022,12,15)]
test_df = df_temp_wind[df_temp_wind.index > datetime(2022,12,15)]


# Define exogenous variables
exog_vars = ['avg_windspd', 'avg_winddir']

# Set the frequency to hourly
train_df= train_df.asfreq('1H') 
#Delete missing values
train_df = train_df.dropna(subset=exog_vars)

# Set Hyper-Parameters
p, d, q = 3, 1, 3
P, D, Q = 0, 0, 0
s = 24

# Fit SARIMAX
sarimax_model = SARIMAX(train_df['avg_temp'],
                order=(p, d, q),
                seasonal_order=(P, D, Q, s),
                exog=train_df[exog_vars])

sarimax_model_fit = sarimax_model.fit()

# Summary
print(sarimax_model_fit.summary())

c:\Users\177119724\AppData\Local\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.

c:\Users\177119724\AppData\Local\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning:

A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.

                               SARIMAX Results                                
==============================================================================
Dep. Variable:               avg_temp   No. Observations:                65024
Model:               SARIMAX(3, 1, 3)   Log Likelihood              157588.046
Date:                Fri, 29 Dec 2023   AIC                        -315158.091
Time:                        18:08:09   BIC                        -315076.349
Sample:                             0   HQIC                       -315132.795
                              - 65024                                         
Covariance Type:                  opg                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
avg_windspd     0.1454      0.004     36.921      0.000       0.138       0.153
avg_winddir -3.687e-05   9.41e-06     -3.918      0.000   -5.53e-05   -1.84e-05
ar.L1           2.4614      0.013    196.384      0.000       2.437       2.486
ar.L2          -1.9693      0.023    -84.233      0.000      -2.015      -1.924
ar.L3           0.5050      0.011     45.665      0.000       0.483       0.527
ma.L1          -0.8106      0.013    -64.192      0.000      -0.835      -0.786
ma.L2          -0.0463      0.004    -10.490      0.000      -0.055      -0.038
ma.L3          -0.0673      0.006    -12.051      0.000      -0.078      -0.056
sigma2          0.0005   1.27e-06    360.363      0.000       0.000       0.000
===================================================================================
Ljung-Box (L1) (Q):                   6.34   Jarque-Bera (JB):            184518.11
Prob(Q):                              0.01   Prob(JB):                         0.00
Heteroskedasticity (H):               1.01   Skew:                            -0.06
Prob(H) (two-sided):                  0.58   Kurtosis:                        11.25
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).


# Get the forecast for the test period
forecast_steps = len(test_df)
forecast_index = test_df.index  # Use the time index of the test dataframe

forecast = sarimax_model_fit.get_forecast(steps=forecast_steps, exog=test_df[exog_vars])

# Plot the actual values
plt.plot(test_df['avg_temp'], label='Actual')

# Plot the forecasted values
plt.plot(forecast_index, forecast.predicted_mean, label='Forecast', color='red')

# Plot confidence intervals
#plt.fill_between(forecast_index, forecast.conf_int()['lower avg_temp'], forecast.conf_int()['upper avg_temp'], color='pink', alpha=0.3)

# Set labels and title
plt.xlabel('Time')
plt.ylabel('Average Temperature')
plt.title('SARIMAX Forecast on Test Data')
plt.legend()

# Show the plot
plt.show()
plt.savefig('../output/Linear Regression Analysis Multiple Predictors SARIMAX Forecast Results.png')

c:\Users\177119724\AppData\Local\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:836: ValueWarning:

No supported index is available. Prediction results will be given with an integer index beginning at `start`.

<Figure size 2000x1500 with 0 Axes>


# Get evaluation data
sarimax_aic = sarimax_model_fit.aic
sarimax_bic = sarimax_model_fit.bic
sarimax_mean_squared_error = sarimax_model_fit.mse
sarimax_sum_squared_error=sarimax_model_fit.sse

print(f'Akaike information criterion | AIC: {sarimax_aic}')
print(f'Bayesian information criterion | BIC: {sarimax_bic}')
print(f'Mean Squared Error | MSE: {sarimax_mean_squared_error}')
print(f'Sum Squared Error | SSE: {sarimax_sum_squared_error}')

Akaike information criterion | AIC: -315158.0914554971
Bayesian information criterion | BIC: -315076.3489885051
Mean Squared Error | MSE: 1.2663238330051534
Sum Squared Error | SSE: 82341.4409173271


# Create a meaningful name with a timestamp
model_name = f"sarimax_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name, 'wb') as f:
    pickle.dump(sarimax_model_fit, f)

print(f"The SARIMAX model has been saved to '{folder_path + model_name}'.")

The SARIMAX model has been saved to '../models/sarimax_model_20231229_1808.pkl'.


#from lazypredict.Supervised import  LazyRegressor
#from sklearn.model_selection import train_test_split

## Define predictor variables
#X = df_temp_wind[['avg_windspd', 'avg_winddir','avg_windgust']]
## Define response variable
#y = df_temp_wind['avg_temp']

#X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2, random_state=123)

#reg = LazyRegressor(verbose=0, ignore_warnings=False, custom_metric=None)
#models,pred = reg.fit(X_train, X_test, y_train, y_test)

#print(models)

## takes too long to run


# Extract features and target
features = ['avg_windspd', 'avg_winddir']
target = 'avg_temp'

X = df_temp_wind[features]
y = df_temp_wind[target]

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Create time-related features
def create_time_feature(data):
        data['dayofmonth'] = data.index.day
        data['dayofweek'] = data.index.dayofweek
        data['quarter'] = data.index.quarter
        data['month'] = data.index.month
        data['year'] = data.index.year
        data['dayofyear'] = data.index.dayofyear
        data['weekofyear'] = data.index.isocalendar().week

create_time_feature(df_temp_wind)
df_temp_wind.head()


from sklearn.model_selection import TimeSeriesSplit
import xgboost as xgb


model = xgb.XGBRegressor(objective='reg:squarederror', n_estimators=1000, max_depth=3, learning_rate=0.1, alpha=0.5)


# Initialize time aware cross-validation
tscv = TimeSeriesSplit(n_splits=5)

# Perform cross-validation
for train_index, test_index in tscv.split(X):
    X_train, X_test = X.iloc[train_index], X.iloc[test_index]
    y_train, y_test = y.iloc[train_index], y.iloc[test_index]
    
    model.fit(X_train, y_train)


test_predictions_index = X_test.index
predictions = model.predict(X_test)

predictions_df = pd.DataFrame(predictions, index=test_predictions_index, columns=['vg_temp_prediction'])


comparison_df = pd.DataFrame({'Actual': y_test, 'Predicted': predictions_df['vg_temp_prediction']})
# Calculate Mean Squared Error (MSE)
mse = mean_squared_error(comparison_df['Actual'], comparison_df['Predicted'])
print(f'Mean Squared Error: {mse}')

# Calculate Mean Absolute Error (MAE)
mae = mean_absolute_error(comparison_df['Actual'], comparison_df['Predicted'])
print(f'Mean Absolute Error: {mae}')

Mean Squared Error: 125.31079889302593
Mean Absolute Error: 9.367344482232422


# 'comparison_df' enthält die tatsächlichen und vorhergesagten Werte
plt.figure(figsize=(12, 6))
plt.plot(comparison_df.index, comparison_df['Actual'], label='Actual',)
plt.plot(comparison_df.index, comparison_df['Predicted'], label='Predicted', )
plt.xlabel('Zeit')
plt.ylabel('vg_temp')
plt.title('Predictions vs. Actual Values')
plt.legend()
plt.show()


# Create a meaningful name with a timestamp
model_name = f"xgb_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name, 'wb') as f:
    pickle.dump(model, f)

print(f"The SARIMAX model has been saved to '{folder_path + model_name}'.")

The SARIMAX model has been saved to '../models/xgb_model_20231229_1808.pkl'.


# Cloning DataFrame to manipulate it for the forecast prediction
forecast_data = pd.DataFrame(df)

# Convert timestamp to datetime format
forecast_data['timestamp'] = pd.to_datetime(forecast_data['run_datetime'])

# Extract day of the year and hour of the day
forecast_data['day_of_year'] = forecast_data['timestamp'].dt.dayofyear
forecast_data['hour'] = forecast_data['timestamp'].dt.hour

# Assuming 'avg_temp' is the target variable
X = forecast_data[['day_of_year', 'hour']]
y = forecast_data['avg_temp']


# Create models
linear_regression_model = make_pipeline(LinearRegression())
gradient_boosting_model = GradientBoostingRegressor()
sgd_regressor_model = make_pipeline(StandardScaler(), SGDRegressor())
svr_model = SVR()


# Split the data into training and testing sets
train_mask = (forecast_data['timestamp'].dt.year >= 2015) & (forecast_data['timestamp'].dt.year <= 2021)
test_mask = (forecast_data['timestamp'].dt.year == 2022)

X_train, X_test = X[train_mask], X[test_mask]
y_train, y_test = y[train_mask], y[test_mask]


# Train models
linear_regression_model_fit = linear_regression_model.fit(X_train, y_train)
gradient_boosting_model_fit = gradient_boosting_model.fit(X_train, y_train)
sgd_regressor_model_fit = sgd_regressor_model.fit(X_train, y_train)
svr_model_fit = svr_model.fit(X_train, y_train)


# Make predictions on the test set
y_pred_linear_regression = linear_regression_model_fit.predict(X_test)
y_pred_gradient_boosting = gradient_boosting_model_fit.predict(X_test)
y_pred_sgd_regressor = sgd_regressor_model_fit.predict(X_test)
y_pred_svr = svr_model_fit.predict(X_test)

# Calculate residuals
residuals_linear_regression = y_test - y_pred_linear_regression
residuals_gradient_boosting = y_test - y_pred_gradient_boosting
residuals_sgd_regressor = y_test - y_pred_sgd_regressor
residuals_svr = y_test - y_pred_svr


# Calculate MSE, MAE, and Mean Residual Error for each model
mse_linear_regression = mean_squared_error(y_test, y_pred_linear_regression)
mae_linear_regression = mean_absolute_error(y_test, y_pred_linear_regression)
mean_residual_error_linear_regression = np.mean(residuals_linear_regression)

mse_gradient_boosting = mean_squared_error(y_test, y_pred_gradient_boosting)
mae_gradient_boosting = mean_absolute_error(y_test, y_pred_gradient_boosting)
mean_residual_error_gradient_boosting = np.mean(residuals_gradient_boosting)

mse_sgd_regressor = mean_squared_error(y_test, y_pred_sgd_regressor)
mae_sgd_regressor = mean_absolute_error(y_test, y_pred_sgd_regressor)
mean_residual_error_sgd_regressor = np.mean(residuals_sgd_regressor)

mse_svr = mean_squared_error(y_test, y_pred_svr)
mae_svr = mean_absolute_error(y_test, y_pred_svr)
mean_residual_error_svr = np.mean(residuals_svr)

# Visualize the predictions, residuals, and tolerance zone
plt.figure(figsize=(15, 20))

# Linear Regression subplot
plt.subplot(4, 2, 1)
plt.plot(forecast_data[test_mask]['timestamp'], y_test, label='Actual')
plt.plot(forecast_data[test_mask]['timestamp'], y_pred_linear_regression, label=f'Linear Regression Predicted\nMSE: {mse_linear_regression:.2f}, MAE: {mae_linear_regression:.2f}, Mean Residual: {mean_residual_error_linear_regression:.2f}')
plt.xlabel('Timestamp')
plt.ylabel('Temperature')
plt.title('Linear Regression')
plt.legend()

# Linear Regression Residuals subplot
plt.subplot(4, 2, 2)
plt.scatter(y_pred_linear_regression, residuals_linear_regression)
plt.axhline(y=0, color='black', linestyle='--', label='Zero Residuals Line')
plt.xlabel('Predicted Temperature')
plt.ylabel('Residuals')
plt.title('Linear Regression Residuals')
plt.legend()

# Gradient Boosting Regressor subplot
plt.subplot(4, 2, 3)
plt.plot(forecast_data[test_mask]['timestamp'], y_test, label='Actual')
plt.plot(forecast_data[test_mask]['timestamp'], y_pred_gradient_boosting, label=f'Gradient Boosting Predicted\nMSE: {mse_gradient_boosting:.2f}, MAE: {mae_gradient_boosting:.2f}, Mean Residual: {mean_residual_error_gradient_boosting:.2f}')
plt.xlabel('Timestamp')
plt.ylabel('Temperature')
plt.title('Gradient Boosting Regressor')
plt.legend()

# Gradient Boosting Residuals subplot
plt.subplot(4, 2, 4)
plt.scatter(y_pred_gradient_boosting, residuals_gradient_boosting)
plt.axhline(y=0, color='black', linestyle='--', label='Zero Residuals Line')
plt.xlabel('Predicted Temperature')
plt.ylabel('Residuals')
plt.title('Gradient Boosting Residuals')
plt.legend()

# SGD Regressor subplot
plt.subplot(4, 2, 5)
plt.plot(forecast_data[test_mask]['timestamp'], y_test, label='Actual')
plt.plot(forecast_data[test_mask]['timestamp'], y_pred_sgd_regressor, label=f'SGD Regressor Predicted\nMSE: {mse_sgd_regressor:.2f}, MAE: {mae_sgd_regressor:.2f}, Mean Residual: {mean_residual_error_sgd_regressor:.2f}')
plt.xlabel('Timestamp')
plt.ylabel('Temperature')
plt.title('SGD Regressor')
plt.legend()

# SGD Regressor Residuals subplot
plt.subplot(4, 2, 6)
plt.scatter(y_pred_sgd_regressor, residuals_sgd_regressor)
plt.axhline(y=0, color='black', linestyle='--', label='Zero Residuals Line')
plt.xlabel('Predicted Temperature')
plt.ylabel('Residuals')
plt.title('SGD Regressor Residuals')
plt.legend()

# SVR subplot
plt.subplot(4, 2, 7)
plt.plot(forecast_data[test_mask]['timestamp'], y_test, label='Actual')
plt.plot(forecast_data[test_mask]['timestamp'], y_pred_svr, label=f'SVR Predicted\nMSE: {mse_svr:.2f}, MAE: {mae_svr:.2f}, Mean Residual: {mean_residual_error_svr:.2f}')
plt.xlabel('Timestamp')
plt.ylabel('Temperature')
plt.title('Support Vector Regressor')
plt.legend()

# SVR Residuals subplot
plt.subplot(4, 2, 8)
plt.scatter(y_pred_svr, residuals_svr)
plt.axhline(y=0, color='black', linestyle='--', label='Zero Residuals Line')
plt.xlabel('Predicted Temperature')
plt.ylabel('Residuals')
plt.title('Support Vector Regressor Residuals')
plt.legend()

plt.tight_layout()
plt.show()
plt.savefig('../output/Linear Regression Analysis Prediction Forecast Results.png')

<Figure size 2000x1500 with 0 Axes>


# Create a meaningful name with a timestamp
model_name_lr_3 = f"lr_model_fit_3_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_gb_3 = f"gb_model_fit_3_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_sgd_3 = f"sgd_model_fit_3_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_svr_3 = f"svr_model_fit_3_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name_lr_3, 'wb') as f:
    pickle.dump(linear_regression_model_fit, f)
with open(folder_path + model_name_gb_3, 'wb') as f:
    pickle.dump(gradient_boosting_model_fit, f)
with open(folder_path + model_name_sgd_3, 'wb') as f:
    pickle.dump(sgd_regressor_model_fit, f)
with open(folder_path + model_name_svr_3, 'wb') as f:
    pickle.dump(svr_model_fit, f)


label1_counts = df['label1'].value_counts()
print(label1_counts)

label1
1    52633
0    12712
Name: count, dtype: int64


# Create Logistic Regressuin Model
logreg_model = LogisticRegression()


# Select the features (independent variable)
X = df[['avg_temp']]

# Select the dependent variable
y = df['label1']

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Standardize the features (reshape to 2D array)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train.values.reshape(-1, 1))
X_test_scaled = scaler.transform(X_test.values.reshape(-1, 1))

# Print the shapes of the sets
print("X_train shape:", X_train_scaled.shape)
print("X_test shape:", X_test_scaled.shape)
print("y_train shape:", y_train.shape)
print("y_test shape:", y_test.shape)

X_train shape: (52276, 1)
X_test shape: (13069, 1)
y_train shape: (52276,)
y_test shape: (13069,)


# Fit Logistic Regression model
logreg_model.fit(X_train_scaled, y_train)

# Add a constant to the features for statsmodels
X_train_sm = sm.add_constant(X_train_scaled)

# Create a logistic regression model using statsmodels
logit_model_sm = sm.Logit(y_train, X_train_sm)

# Fit the logit sm  model
logit_model = logit_model_sm.fit()

Optimization terminated successfully.
         Current function value: 0.375341
         Iterations 7


from yellowbrick.classifier import ConfusionMatrix, ClassificationReport

# Print AIC
aic_value = logit_model.aic
print("AIC:", aic_value)
print("Label 0: Extreme Weather Event \n Label 1: Blue Sky Day")

# Make predictions on the test set using the result object
X_test_sm = sm.add_constant(X_test_scaled)
y_prob = logit_model.predict(X_test_sm)  # Probabilities
y_pred = (y_prob > 0.5).astype(int)  # Convert probabilities to binary predictions

# Confusion Matrix as Diagram with custom colors
cm_viz = ConfusionMatrix(logreg_model, classes=[0, 1], cmap='RdYlGn')
cm_viz.score(X_test_scaled, y_test)
cm_viz.show()
plt.show()
plt.savefig('../output/Logistic Regression Predicting WEP by Temperature Confussion Matrix.png')

print("\nClassification Report:")
print(classification_report(y_test, y_pred))

AIC: 39246.639678249965
Label 0: Extreme Weather Event 
 Label 1: Blue Sky Day

Classification Report:
              precision    recall  f1-score   support

           0       0.39      0.20      0.26      2515
           1       0.83      0.93      0.88     10554

    accuracy                           0.79     13069
   macro avg       0.61      0.56      0.57     13069
weighted avg       0.75      0.79      0.76     13069

<Figure size 800x550 with 0 Axes>


# Visualize results
plt.figure(figsize=(12, 6))

# Plot the data points on the logistic regression curve based on probabilities
y_vals_prob = logit_model.predict(sm.add_constant(scaler.transform(X_test.values.reshape(-1, 1))))
plt.scatter(X_test, y_test, color='blue', label='Actual Labels, yi')
plt.scatter(X_test, y_vals_prob, color='red', marker='o', label='Predicted Probabilities, p(Xi)')

# Plot the linear part of the logistic regression curve
x_vals = np.linspace(X_test.min(), X_test.max(), 100)
linear_combination = logit_model.params[0] + logit_model.params[1] * scaler.transform(x_vals.reshape(-1, 1))

# Logistic function (sigmoid) to convert linear combination to probability
logistic_function = 1 / (1 + np.exp(-linear_combination))

plt.plot(x_vals, linear_combination, color='green', linewidth=2, label='Logit, f(x)')
plt.plot(x_vals, logistic_function, color='orange', linewidth=2, label='Logistic Regression Curve, p(x)')

# Plot the threshold line at y=0.5
plt.axhline(y=0.5, color='purple', linestyle='--', label='Threshold = 0.5')

# Set y-axis limits
plt.ylim(-0.01, 1.01)

# Display AIC value
plt.text(plt.xlim()[1], plt.ylim()[1], f'AIC = {aic_value:.2f}', ha='right', va='top', fontsize=10, color='purple')

plt.xlabel('Temperature')
plt.ylabel('Probability of Extreme Event')
plt.suptitle('Logistic Regression Analysis')
plt.title('Label 0: Extreme Weather Event \n Label 1: Blue Sky Day')
plt.legend()
plt.show()
plt.savefig('../output/Logistic Regression Analysis Predicting WEP by Temperature Results.png')

<Figure size 800x550 with 0 Axes>


# Create a meaningful name with a timestamp
model_name_logreg= f"logreg_model_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name_logreg, 'wb') as f:
    pickle.dump(logreg_model, f)


# show data as pie chart
labels = ['Blue Sky', 'Extreme Weather']
sizes = df['label1'].value_counts()
colors = [ 'lightskyblue','lightcoral']
explode = (0.1, 0)  # explode first


# Plotting the pie chart
plt.figure(figsize=(6, 6))
plt.pie(sizes, explode=explode, labels=labels, colors=colors, autopct=lambda p: '{:.0f} ({:.1f}%)'.format(p * sum(sizes) / 100, p),
        shadow=True, startangle=140)
plt.axis('equal')
plt.title('Blue Sky Events vs. Extreme Weather Events')
plt.show()
plt.savefig('../output/Binary Classification Blue Sky vs Extreme Weather Events.png')

<Figure size 800x550 with 0 Axes>


#copy df
df_class =df

# Calculate the mean value of the column
mean_value = df_class["avg_pressure_change"].mean()

#set "avg_temp_change" to 0 in the first row, as there is no temp change
df_class.at[0, 'avg_temp_change'] = 0

#drop columns, which are not numeric and not needed
df_class = df_class.drop(columns = ['wep','wind_direction_label', 'run_datetime','label2','avg_pressure_change','label0', 'display_label'])

#show adjusted df
df_class.head()


df_class["label1"]

0        1
1        1
2        1
3        1
4        1
        ..
65340    0
65341    1
65342    1
65343    1
65344    1
Name: label1, Length: 65345, dtype: int64


# checking for missing values in the data
df_class.isnull().sum()

avg_temp                 0
avg_temp_celsius         0
min_wet_bulb_temp        0
avg_dewpoint             0
avg_temp_change          0
avg_windspd              0
max_windgust             0
avg_winddir              0
avg_winddir_sin          0
avg_winddir_cos          0
max_cumulative_precip    0
max_snow_density_6       0
max_cumulative_snow      0
max_cumulative_ice       0
label1                   0
dtype: int64


#from lazypredict.Supervised import LazyClassifier

## Define predictor variables
X = df_class.drop(columns = ['label1', 'max_snow_density_6'])

## Define response variable
y = df_class['label1']

#X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

#clf = LazyClassifier(verbose=0, ignore_warnings=False, custom_metric=None)
#models, pred = clf.fit(X_train, X_test, y_train, y_test)

#print(models)

## take too long to run


# Split the dataset into training and testing data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Initialize the classifiers
classifiers = {
    'ExtraTrees': ExtraTreesClassifier(),
    'XGBoost': XGBClassifier(),
    'LGBM': LGBMClassifier(verbose=-1),
    'RandomForest': RandomForestClassifier()
}


# Train the classifiers and collect accuracies
accuracies = {}

for idx, (name, clf) in enumerate(classifiers.items()):
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    accuracy = accuracy_score(y_test, y_pred)
    accuracies[name] = accuracy
    print(f'{name} Accuracy: {accuracy}')

ExtraTrees Accuracy: 0.9931134746346316
XGBoost Accuracy: 0.9941847119136888
LGBM Accuracy: 0.9935725763256561
RandomForest Accuracy: 0.9935725763256561


# Set up subplots
fig, axes = plt.subplots(1, len(classifiers), figsize=(18, 5))

# Iterate over classifiers and plot confusion matrices
for ax, (name, model) in zip(axes, classifiers.items()):
    # Create a pipeline with StandardScaler and the current model
    pipe = make_pipeline(StandardScaler(), model)
    pipe.fit(X_train, y_train)
    y_pred = pipe.predict(X_test)

    # Calculate confusion matrix
    cm = confusion_matrix(y_test, y_pred)

    # Plot heatmap
    sns.heatmap(cm, annot=True, fmt=".0f", linewidths=.5, square=True, cmap='Blues', ax=ax)

    # Set labels and title
    ax.set_xlabel('Predicted Label')
    ax.set_ylabel('True Label')
    ax.set_title(f'{name} Confusion Matrix')

# Adjust layout for better spacing
plt.tight_layout()

# Show the plot
plt.show()
plt.savefig('../output/Binary Classification Confussion Matrix.png')

<Figure size 800x550 with 0 Axes>


# Iterate over classifiers and plot confusion matrices
for ax, (name, model) in zip(axes, classifiers.items()):
    # Create a pipeline with StandardScaler and the current model
    pipe = make_pipeline(StandardScaler(), model)
    pipe.fit(X_train, y_train)
    y_pred = pipe.predict(X_test)

    # Calculate classification report with more decimal places
    report = classification_report(y_test, y_pred, digits=4)

    print(name, ":\n", report)

ExtraTrees :
               precision    recall  f1-score   support

           0     0.9911    0.9742    0.9826      2515
           1     0.9939    0.9979    0.9959     10554

    accuracy                         0.9933     13069
   macro avg     0.9925    0.9860    0.9892     13069
weighted avg     0.9933    0.9933    0.9933     13069

XGBoost :
               precision    recall  f1-score   support

           0     0.9931    0.9753    0.9842      2515
           1     0.9942    0.9984    0.9963     10554

    accuracy                         0.9940     13069
   macro avg     0.9936    0.9869    0.9902     13069
weighted avg     0.9940    0.9940    0.9939     13069

LGBM :
               precision    recall  f1-score   support

           0     0.9947    0.9726    0.9835      2515
           1     0.9935    0.9988    0.9961     10554

    accuracy                         0.9937     13069
   macro avg     0.9941    0.9857    0.9898     13069
weighted avg     0.9937    0.9937    0.9937     13069

RandomForest :
               precision    recall  f1-score   support

           0     0.9947    0.9714    0.9829      2515
           1     0.9932    0.9988    0.9960     10554

    accuracy                         0.9935     13069
   macro avg     0.9940    0.9851    0.9894     13069
weighted avg     0.9935    0.9935    0.9935     13069


for ax, (name, model) in zip(axes, classifiers.items()):

    model_name= f"{name}_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

    # Specify the folder path
    folder_path = "../models/"

    # Save the model to the specified folder with the created name
    with open(folder_path + model_name, 'wb') as f:
        pickle.dump(model, f)


# Creating new dataframe that only includes extreme weather observations
df_extreme_weather = pd.DataFrame(df)
df_extreme_weather = df_extreme_weather[df_extreme_weather['label1'] != 1]

# Dropping NaN
df_extreme_weather = df_extreme_weather.dropna()
df_extreme_weather


knn_model = KNeighborsClassifier(
    n_neighbors=5, 
    weights='uniform', 
    algorithm='auto', 
    leaf_size=30, 
    p=2, 
    metric='minkowski', 
    metric_params=None, 
    n_jobs=None)

svm_model = svm.SVC(
    C=1.0, 
    kernel='rbf', 
    degree=3, 
    gamma='scale', 
    coef0=0.0, 
    tol=0.001, 
    cache_size=800, 
    class_weight='balanced', 
    random_state=42)

dtc_model = DecisionTreeClassifier(
    criterion="entropy", 
    class_weight="balanced", 
    max_depth=10, 
    min_samples_leaf=3, 
    max_features=None, 
    random_state=42)

gbc_model = GradientBoostingClassifier(
    loss='log_loss', 
    learning_rate=0.1, 
    n_estimators=100, 
    subsample=1.0, 
    criterion='friedman_mse', 
    min_samples_split=2, 
    min_samples_leaf=3, 
    min_weight_fraction_leaf=0.0, 
    max_depth=10, 
    )


# Visual inspection of the class distribution
wep_counts = df_extreme_weather['wep'].value_counts()
print(wep_counts)
sns.countplot(x = 'wep', data=df_extreme_weather, palette = "Set1")

wep
Mild snowfall                                          3892
Moderate snowfall                                      2483
Moderate rain                                          2322
Heavy snowfall with accumulated snow                    926
Frondurchlauf / Continuous freezing rain                875
Storm with freezing rain / Heavy snow- and icestorm     401
Snowstorm with high precipitation                       210
Name: count, dtype: int64

<Axes: xlabel='wep', ylabel='count'>


# Defining the dependent and independent variables
X = df_extreme_weather.drop(columns=["label0", "label1", "label2", "wep", "avg_winddir", "wind_direction_label", "avg_temp_celsius", "run_datetime","display_label"])# Features
y = df_extreme_weather['label2'] # Labels


# Split into train and test
X_train, X_test, Y_train, Y_test = train_test_split(X, y, test_size=0.3, random_state=42)

print('X_train count: ', X_train.shape[0])
print('X_test count: ', X_test.shape[0])
print('Y_train count: ', Y_train.shape[0])
print('Y_test count: ', Y_test.shape[0])

X_train count:  7776
X_test count:  3333
Y_train count:  7776
Y_test count:  3333


knn_model_fit = knn_model.fit(X_train, Y_train)
svm_model_fit = svm_model.fit(X_train, Y_train)
dtc_model_fit = dtc_model.fit(X_train, Y_train)
gbc_model_fit = gbc_model.fit(X_train, Y_train)


# Predict using the trained KNN model on the test set
Y_pred_knn = knn_model_fit.predict(X_test)

# Generate the confusion matrix
cm_knn = confusion_matrix(Y_test, Y_pred_knn)

# Visualize the Confusion Matrix
plt.figure(figsize=(10, 7))
sns.heatmap(cm_knn, annot=True, fmt='d', cmap='Blues')
plt.title('K-Nearest Neighbors Confusion Matrix ')
plt.xlabel('Predicted')
plt.ylabel('True')
plt.show()
plt.savefig('../output/Multiclass Classification KNN Confusion Matrix.png')

# Display the classification report
print("Classification Report:\n",classification_report(Y_test, Y_pred_knn))

              precision    recall  f1-score   support

         0.0       0.88      0.88      0.88       741
         1.0       0.78      0.82      0.80       274
         2.0       0.88      0.78      0.82       120
         3.0       1.00      1.00      1.00       679
         4.0       0.97      0.85      0.90        66
         5.0       0.95      0.97      0.96      1190
         6.0       0.84      0.80      0.82       263

    accuracy                           0.92      3333
   macro avg       0.90      0.87      0.88      3333
weighted avg       0.92      0.92      0.92      3333

<Figure size 800x550 with 0 Axes>


# Predict using the trained SVM model on the test set
Y_pred_svm = svm_model_fit.predict(X_test)

# Generate the confusion matrix
cm_svm = confusion_matrix(Y_test, Y_pred_svm)

# Visualize the Confusion Matrix
plt.figure(figsize=(10, 7))
sns.heatmap(cm_svm, annot=True, fmt='d', cmap='Blues')
plt.title('Support Vector Machine Confusion Matrix')
plt.xlabel('Predicted')
plt.ylabel('True')
plt.show()
plt.savefig('../output/Multiclass Classification SVM Confusion Matrix.png')

# Display the classification report
print("Classification Report:\n",classification_report(Y_test, Y_pred_svm))

              precision    recall  f1-score   support

         0.0       0.81      0.61      0.69       741
         1.0       0.41      0.69      0.52       274
         2.0       0.78      0.55      0.64       120
         3.0       0.95      1.00      0.98       679
         4.0       0.70      0.94      0.80        66
         5.0       0.94      0.92      0.93      1190
         6.0       0.68      0.67      0.68       263

    accuracy                           0.82      3333
   macro avg       0.75      0.77      0.75      3333
weighted avg       0.84      0.82      0.82      3333

<Figure size 800x550 with 0 Axes>


from sklearn.tree import plot_tree

# Make predictions using the Decision Tree Classifier (DTC) model on the test set
Y_pred_dtc = dtc_model_fit.predict(X_test)

# Calculate the confusion matrix
cm_dtc = confusion_matrix(Y_test, Y_pred_dtc)

# Visualize the Confusion Matrix using seaborn heatmap
plt.figure(figsize=(10,7))
sns.heatmap(cm_dtc, annot=True, fmt='d', cmap='Blues')  
plt.title('Decision Tree Classifier confusion matrix ')
plt.xlabel('Predicted')
plt.ylabel('True')
plt.show()

# Save the confusion matrix plot
plt.savefig('../output/Multiclass Classification DTC Confussion Matrix.png')

# Print the classification report
print("Classification Report:\n",classification_report(Y_test, Y_pred_dtc))

              precision    recall  f1-score   support

         0.0       0.94      0.92      0.93       741
         1.0       0.84      0.93      0.88       274
         2.0       0.89      0.91      0.90       120
         3.0       1.00      0.99      1.00       679
         4.0       0.95      0.94      0.95        66
         5.0       0.98      0.97      0.98      1190
         6.0       0.92      0.92      0.92       263

    accuracy                           0.95      3333
   macro avg       0.93      0.94      0.94      3333
weighted avg       0.96      0.95      0.95      3333

<Figure size 800x550 with 0 Axes>


# Visualize the Decision Tree using plot_tree from sklearn
plt.figure(figsize=(40,20))
fig = plot_tree(dtc_model,
                filled = True,
                rounded = True,
                fontsize = 14)


# Make predictions using the Gradient Boosting Classifier (GBC) model on the test set
Y_pred_gbc = gbc_model_fit.predict(X_test)

# Calculate the confusion matrix
cm_gbc = confusion_matrix(Y_test, Y_pred_gbc)

# Visualize the Confusion Matrix using seaborn heatmap
plt.figure(figsize=(10, 7))
sns.heatmap(cm_gbc, annot=True, fmt='d', cmap='Blues')  
plt.title('Gradient Boosting Classifier confusion matrix ')
plt.xlabel('Predicted')
plt.ylabel('True')
plt.show()

# Save the confusion matrix plot
plt.savefig('../output/Multiclass Classification GBC Confusion Matrix.png')

# Print the classification report
print("Classification Report:\n",classification_report(Y_test, Y_pred_gbc))

# Visualize the Decision Tree using plot_tree from sklearn (assuming it's a GBC model)
plt.figure(figsize=(16, 9))

              precision    recall  f1-score   support

         0.0       0.97      0.97      0.97       741
         1.0       0.93      0.96      0.94       274
         2.0       0.96      0.92      0.94       120
         3.0       1.00      1.00      1.00       679
         4.0       0.98      0.94      0.96        66
         5.0       0.99      0.99      0.99      1190
         6.0       0.95      0.95      0.95       263

    accuracy                           0.98      3333
   macro avg       0.97      0.96      0.96      3333
weighted avg       0.98      0.98      0.98      3333

<Figure size 26500x11500 with 0 Axes>

<Figure size 800x550 with 0 Axes>

<Figure size 26500x11500 with 0 Axes>


# Create a meaningful name with a timestamp
model_name_knn = f"knn_model_fit_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_svm = f"svm_model_fit_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_dtc = f"dtc_model_fit_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"
model_name_gbc = f"gbc_model_fit_{datetime.now().strftime('%Y%m%d_%H%M')}.pkl"

# Specify the folder path
folder_path = "../models/"

# Save the model to the specified folder with the created name
with open(folder_path + model_name_knn, 'wb') as f:
    pickle.dump(knn_model_fit, f)
with open(folder_path + model_name_svm, 'wb') as f:
    pickle.dump(svm_model_fit, f)
with open(folder_path + model_name_dtc, 'wb') as f:
    pickle.dump(dtc_model_fit, f)
with open(folder_path + model_name_gbc, 'wb') as f:
    pickle.dump(gbc_model_fit, f)


!pip freeze requirements.txt

aiohttp==3.9.1
aiosignal==1.3.1
alembic==1.13.0
altair==5.2.0
annotated-types==0.6.0
anyio==3.7.1
asttokens==2.4.1
attrs==23.1.0
certifi==2023.11.17
charset-normalizer==3.3.2
click==8.1.7
colorama==0.4.6
colorlog==6.8.0
comm==0.2.0
contourpy==1.2.0
cycler==0.12.1
Cython==3.0.6
dataclasses-json==0.6.3
debugpy==1.8.0
decorator==5.1.1
distlib==0.3.7
executing==2.0.1
fastapi==0.104.1
fastjsonschema==2.19.0
filelock==3.13.1
fonttools==4.45.1
frozenlist==1.4.0
greenlet==3.0.1
idna==3.6
imbalanced-learn==0.11.0
imblearn==0.0
ipykernel==6.26.0
ipython==8.17.2
jedi==0.19.1
Jinja2==3.1.2
joblib==1.3.2
jsonpatch==1.33
jsonpointer==2.4
jsonschema==4.20.0
jsonschema-specifications==2023.11.1
jupyter_client==8.6.0
jupyter_core==5.5.0
kiwisolver==1.4.5
langchain==0.0.346
langchain-core==0.0.10
langsmith==0.0.69
lazypredict==0.2.12
lightgbm==4.1.0
Mako==1.3.0
MarkupSafe==2.1.3
marshmallow==3.20.1
matplotlib==3.8.2
matplotlib-inline==0.1.6
multidict==6.0.4
mypy-extensions==1.0.0
nbformat==5.9.2
nest-asyncio==1.5.8
numpy==1.26.2
opensearch-py==2.4.2
optuna==3.4.0
packaging==23.2
pandas==2.1.3
parso==0.8.3
patsy==0.5.4
Pillow==10.1.0
platformdirs==4.0.0
plotly==5.18.0
pmdarima==2.0.4
podman==4.8.0.post1
podman-compose==1.0.6
prompt-toolkit==3.0.40
protobuf==4.25.1
psutil==5.9.6
pure-eval==0.2.2
pyarrow==14.0.1
pydantic==2.5.2
pydantic_core==2.14.5
Pygments==2.16.1
pyparsing==3.1.1
PyPDF2==3.0.1
python-dateutil==2.8.2
python-dotenv==1.0.0
pytz==2023.3.post1
pywin32==306
pyxdg==0.28
PyYAML==6.0.1
pyzmq==25.1.1
referencing==0.31.0
requests==2.31.0
rpds-py==0.13.1
scikit-learn==1.3.2
scipy==1.11.4
seaborn==0.13.0
setuptools==69.0.2
six==1.16.0
skforecast==0.11.0
sniffio==1.3.0
SQLAlchemy==2.0.23
stack-data==0.6.3
starlette==0.27.0
statsmodels==0.14.0
tenacity==8.2.3
threadpoolctl==3.2.0
toolz==0.12.0
tornado==6.3.3
tqdm==4.66.1
traitlets==5.13.0
typing-inspect==0.9.0
typing_extensions==4.8.0
tzdata==2023.3
urllib3==2.1.0
vegafusion==1.4.5
vegafusion-python-embed==1.4.5
virtualenv==20.25.0
vl-convert-python==1.1.0
wcwidth==0.2.10
xgboost==2.0.2
yarl==1.9.4
yellowbrick==1.5

	count	mean	min	25%	50%	75%	max	std
run_datetime	65345	2019-04-06 14:09:11.362766848	2015-07-15 00:00:00	2017-05-25 23:00:00	2019-04-07 01:00:00	2021-02-14 16:00:00	2022-12-27 08:00:00	NaN
avg_temp	65345.0	279.574328	243.849393	271.114219	279.882735	289.903226	300.934144	11.383325
avg_temp_celsius	65345.0	6.424328	-29.300607	-2.035781	6.732735	16.753226	27.784144	11.383325
min_wet_bulb_temp	65345.0	273.601253	239.056465	265.232651	274.245761	284.128378	294.642287	11.914203
avg_dewpoint	65345.0	274.972135	238.429236	267.161062	275.429991	284.366793	294.97882	10.93508
avg_temp_change	65344.0	-0.000361	-1.168328	-0.095711	0.012257	0.10334	1.104699	0.191027
avg_windspd	65345.0	2.440013	0.643393	1.882336	2.34648	2.893087	6.470306	0.77288
max_windgust	65345.0	12.659095	3.918856	9.966063	12.266102	14.895518	28.236084	3.721124
avg_winddir	65345.0	210.247462	21.17728	165.321173	222.129253	262.922675	346.805959	70.896648
avg_winddir_sin	65345.0	-0.008188	-1.0	-0.715716	-0.012884	0.699482	1.0	0.707159
avg_winddir_cos	65345.0	0.003304	-1.0	-0.704408	0.008696	0.708963	1.0	0.70701
max_cumulative_precip	65345.0	4.796664	0.0	0.1	1.528	6.2	83.471	7.64858
max_snow_density_6	65345.0	21.123859	0.0	0.0	0.0	0.0	722.582771	67.439046
max_cumulative_snow	65345.0	17.454941	0.0	0.0	0.0	8.257	920.464	53.768892
max_cumulative_ice	65345.0	0.056319	0.0	0.0	0.0	0.0	9.955316	0.405855
avg_pressure_change	56672.0	0.017391	-168.315674	-19.263538	-0.071299	19.549097	153.811455	33.88573
label0	65345.0	2.027072	0.0	0.0	2.0	3.0	7.0	1.933712
label1	65345.0	0.805463	0.0	1.0	1.0	1.0	1.0	0.395847
label2	12712.0	3.06191	0.0	1.0	3.0	5.0	6.0	2.126446

	avg_temp	avg_windspd	avg_winddir	avg_winddir_sin	avg_winddir_cos	max_cumulative_precip	label1
0	287.389224	3.386380	80.302464	-0.981653	0.190676	2.009	1
1	287.378997	3.326687	76.866373	0.994736	0.102466	1.209	1
2	287.388845	3.243494	76.258867	0.758262	0.651950	0.400	1
3	287.427324	3.145505	78.299616	0.237897	-0.971290	0.000	1
4	287.489158	3.047607	84.632852	0.189006	-0.981976	0.000	1

	avg_temp	avg_winddir	avg_windspd	avg_winddir_sin	avg_winddir_cos	max_cumulative_precip	max_snow_density_6	max_cumulative_snow	max_cumulative_ice
run_datetime
2015-07-15 00:00:00	0.000000	80.302464	3.386380	-0.981653	0.190676	2.009	0.0	0.0	0.0
2015-07-15 01:00:00	287.378997	76.866373	3.326687	0.994736	0.102466	1.209	0.0	0.0	0.0
2015-07-15 02:00:00	287.388845	76.258867	3.243494	0.758262	0.651950	0.400	0.0	0.0	0.0
2015-07-15 03:00:00	287.427324	78.299616	3.145505	0.237897	-0.971290	0.000	0.0	0.0	0.0
2015-07-15 04:00:00	287.489158	84.632852	3.047607	0.189006	-0.981976	0.000	0.0	0.0	0.0

	avg_temp	avg_windspd	avg_winddir	avg_windgust
run_datetime
2015-07-15 00:00:00	287.389224	3.386380	80.302464	11.136197
2015-07-15 01:00:00	287.378997	3.326687	76.866373	11.002795
2015-07-15 02:00:00	287.388845	3.243494	76.258867	10.700595
2015-07-15 03:00:00	287.427324	3.145505	78.299616	10.323983
2015-07-15 04:00:00	287.489158	3.047607	84.632852	9.921157

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	2.4725	0.021	116.934	0.000	2.431	2.514
ar.L2	-2.0030	0.038	-52.033	0.000	-2.078	-1.928
ar.L3	0.5260	0.018	29.575	0.000	0.491	0.561
ma.L1	-0.8029	0.021	-37.572	0.000	-0.845	-0.761
ma.L2	-0.0242	0.005	-5.296	0.000	-0.033	-0.015
ma.L3	-0.0307	0.006	-5.555	0.000	-0.042	-0.020
sigma2	0.0005	1.2e-06	384.413	0.000	0.000	0.000

	Unnamed: 0	substation	run_datetime	valid_datetime	avg_temp	avg_windspd	avg_windgust	avg_pressure_change	...	MA_avg_winddir_month	MA_avg_temp_change_month	label0	label1	label2	label3	wep	storm_id2	year	month
0	0	Bancroft	2015-07-15 00:00:00	2015-07-15 00:00:00	287.389224	3.386380	11.136197	52.892217	...	80.302464	NaN	0	1	NaN	1	Blue sky day	NaN	2015	7
1	1	Bancroft	2015-07-15 01:00:00	2015-07-15 01:00:00	287.378997	3.326687	11.002795	50.256685	...	78.584418	-0.010227	0	1	NaN	1	Blue sky day	NaN	2015	7
2	2	Bancroft	2015-07-15 02:00:00	2015-07-15 02:00:00	287.388845	3.243494	10.700595	47.944054	...	77.809235	-0.000189	3	1	NaN	1	Blue sky day	NaN	2015	7
3	3	Bancroft	2015-07-15 03:00:00	2015-07-15 03:00:00	287.427324	3.145505	10.323983	45.855264	...	77.931830	0.012700	2	1	NaN	1	Blue sky day	NaN	2015	7
4	4	Bancroft	2015-07-15 04:00:00	2015-07-15 04:00:00	287.489158	3.047607	9.921157	44.823453	...	79.272034	0.024983	2	1	NaN	1	Blue sky day	NaN	2015	7

	Name	Description	Role	Type	Format
0	run_datetime	Date and time when the weather observations we...	ID / predictor	numerical continuous / ID	<class 'pandas._libs.tslibs.timestamps.Timesta...
1	wep	Weather Event Type (WEP) is a categorization o...	response	categorical nominal	<class 'str'>
2	avg_temp	The average temperature measured at two meters...	response / predictor	numerical continuous	<class 'numpy.float64'>
3	min_wet_bulb_temp	Minimum wet bulb temperature recorded during t...	predictor	numerical continuous	<class 'numpy.float64'>
4	avg_dewpoint	Average dewpoint temperature observed during t...	predictor	numerical continuous	<class 'numpy.float64'>
5	avg_temp_change	Average change in temperature during the obser...	predictor	numerical continuous	<class 'numpy.float64'>
6	avg_windspd	Average wind speed measured during the recordi...	predictor	numerical continuous	<class 'numpy.float64'>
7	max_windgust	Maximum wind gust observed during the recordin...	predictor	numerical continuous	<class 'numpy.float64'>
8	avg_winddir	Average wind direction (in degree) observed du...	predictor	numerical continuous	<class 'numpy.float64'>
9	wind_direction_label	Wind direction (in cardinal direction) observe...	predictor	categorical ordinal	<class 'str'>
10	max_cumulative_precip	Maximum cumulative precipitation recorded, con...	predictor	numerical continuous	<class 'numpy.float64'>
11	max_snow_density_6	Maximum snow density at a depth of 6 inches, c...	predictor	numerical continuous	<class 'numpy.float64'>
12	max_cumulative_snow	Maximum cumulative snow recorded, considering ...	predictor	numerical continuous	<class 'numpy.float64'>
13	max_cumulative_ice	Maximum cumulative ice recorded, considering a...	predictor	numerical continuous	<class 'numpy.float64'>
14	avg_pressure_change	Average change in atmospheric pressure during ...	predictor	numerical continuous	<class 'numpy.float64'>

Dep. Variable:	y	No. Observations:	65320
Model:	SARIMAX(3, 1, 3)	Log Likelihood	157921.228
Date:	Fri, 29 Dec 2023	AIC	-315828.456
Time:	18:07:02	BIC	-315764.847
Sample:	0	HQIC	-315808.775
	- 65320
Covariance Type:	opg

Ljung-Box (L1) (Q):	4.99	Jarque-Bera (JB):	230383.67
Prob(Q):	0.03	Prob(JB):	0.00
Heteroskedasticity (H):	1.00	Skew:	-0.09
Prob(H) (two-sided):	0.91	Kurtosis:	12.20

Model	Adjusted R-Squared	R-Squared	RMSE	Time Taken
XGBRegressor	0.13	0.13	10.66	0.44
HistGradientBoostingRegressor	0.12	0.12	10.72	0.80
LGBMRegressor	0.12	0.12	10.72	0.17
RandomForestRegressor	0.12	0.12	10.72	69.52
ExtraTreesRegressor	0.10	0.10	10.83	14.12
MLPRegressor	0.08	0.08	10.95	30.44
GradientBoostingRegressor	0.08	0.08	10.95	14.94
NuSVR	0.07	0.07	11.02	127.68
SVR	0.06	0.06	11.09	129.99
KNeighborsRegressor	0.05	0.05	11.11	0.11
BaggingRegressor	0.05	0.05	11.13	7.32
AdaBoostRegressor	0.03	0.03	11.25	1.34
LassoCV	0.01	0.01	11.34	0.42
LassoLarsCV	0.01	0.01	11.34	0.06
Lars	0.01	0.01	11.34	0.17
LarsCV	0.01	0.01	11.34	0.15
TransformedTargetRegressor	0.01	0.01	11.34	0.03
LassoLarsIC	0.01	0.01	11.34	0.04
LinearRegression	0.01	0.01	11.34	0.02
OrthogonalMatchingPursuitCV	0.01	0.01	11.34	0.05
Ridge	0.01	0.01	11.34	0.02
RidgeCV	0.01	0.01	11.34	0.03
ElasticNetCV	0.01	0.01	11.34	0.30
BayesianRidge	0.01	0.01	11.34	0.03
PoissonRegressor	0.01	0.01	11.34	0.05
SGDRegressor	0.01	0.01	11.34	0.07
OrthogonalMatchingPursuit	0.01	0.01	11.35	0.02
TweedieRegressor	0.01	0.01	11.35	0.02
GammaRegressor	0.01	0.01	11.35	0.14
ElasticNet	0.01	0.01	11.36	0.03
LassoLars	0.01	0.01	11.39	0.02
Lasso	0.01	0.01	11.39	0.04
HuberRegressor	0.00	0.01	11.39	0.14
LinearSVR	0.00	0.00	11.42	0.10
DummyRegressor	-0.00	-0.00	11.42	0.01
PassiveAggressiveRegressor	-0.20	-0.20	12.51	0.04
RANSACRegressor	-0.48	-0.47	13.87	0.13
DecisionTreeRegressor	-0.64	-0.64	14.64	0.86
ExtraTreeRegressor	-0.69	-0.68	14.82	0.15

	avg_temp	avg_temp_celsius	min_wet_bulb_temp	avg_dewpoint	avg_temp_change	avg_windspd	max_windgust	avg_winddir	avg_winddir_sin	avg_winddir_cos	max_cumulative_precip	label1
0	287.389224	14.239224	280.809506	280.735246	0.000000	3.386380	14.899891	80.302464	-0.981653	0.190676	2.009	1
1	287.378997	14.228997	280.809506	280.414058	-0.010227	3.326687	14.899891	76.866373	0.994736	0.102466	1.209	1
2	287.388845	14.238845	280.809506	280.187074	0.009848	3.243494	14.899891	76.258867	0.758262	0.651950	0.400	1
3	287.427324	14.277324	280.809506	280.049330	0.038479	3.145505	14.899891	78.299616	0.237897	-0.971290	0.000	1
4	287.489158	14.339158	280.809506	279.980697	0.061834	3.047607	14.702229	84.632852	0.189006	-0.981976	0.000	1

Model	Accuracy	Balanced Accuracy	ROC AUC	F1 Score	Time Taken
ExtraTreesClassifier	0.99	0.99	0.99	0.99	4.13
XGBClassifier	0.99	0.99	0.99	0.99	0.63
RandomForestClassifier	0.99	0.99	0.99	0.99	29.54
LGBMClassifier	0.99	0.99	0.99	0.99	0.67
BaggingClassifier	0.99	0.98	0.98	0.99	11.56
DecisionTreeClassifier	0.99	0.98	0.98	0.99	1.54
AdaBoostClassifier	0.99	0.98	0.98	0.99	14.65
LinearSVC	0.99	0.97	0.97	0.99	1.62
CalibratedClassifierCV	0.99	0.97	0.97	0.99	1.48
LogisticRegression	0.98	0.97	0.97	0.98	0.43
QuadraticDiscriminantAnalysis	0.98	0.97	0.97	0.98	0.24
GaussianNB	0.97	0.97	0.97	0.97	0.11
SVC	0.99	0.97	0.97	0.99	22.36
SGDClassifier	0.99	0.97	0.97	0.99	0.16
ExtraTreeClassifier	0.97	0.96	0.96	0.97	0.12
Perceptron	0.98	0.95	0.95	0.98	0.20
PassiveAggressiveClassifier	0.98	0.95	0.95	0.98	0.23
KNeighborsClassifier	0.97	0.92	0.92	0.96	7.63
BernoulliNB	0.89	0.92	0.92	0.90	0.14
NearestCentroid	0.83	0.86	0.86	0.85	0.13
LinearDiscriminantAnalysis	0.91	0.76	0.76	0.90	0.37
RidgeClassifier	0.89	0.72	0.72	0.87	0.11
RidgeClassifierCV	0.89	0.72	0.72	0.87	0.13
DummyClassifier	0.81	0.50	0.50	0.72	0.06

	run_datetime	wep	avg_temp	avg_temp_celsius	min_wet_bulb_temp	avg_dewpoint	avg_temp_change	avg_windspd	max_windgust	avg_winddir	...	wind_direction_label	max_cumulative_precip	max_snow_density_6	max_cumulative_snow	max_cumulative_ice	avg_pressure_change	label0	label1	label2	display_label
2238	2015-10-16 13:00:00	Moderate rain	276.935854	3.785854	271.085195	273.479160	-0.196776	2.610431	12.095491	302.159235	...	Northwest	6.115	0.000000	24.979	0.0	25.754899	4	0	3.0	Moderate rain
2239	2015-10-16 14:00:00	Moderate rain	276.763328	3.613328	271.085195	273.237344	-0.172526	2.657227	12.095491	302.695638	...	Northwest	6.215	0.000000	25.897	0.0	26.257896	4	0	3.0	Moderate rain
2240	2015-10-16 15:00:00	Moderate rain	276.597296	3.447296	271.085195	272.988734	-0.166032	2.700702	12.095491	303.377860	...	Northwest	6.315	0.000000	26.591	0.0	25.630187	4	0	3.0	Moderate rain
2241	2015-10-16 16:00:00	Moderate rain	276.430385	3.280385	271.085195	272.769039	-0.166911	2.743809	12.095491	304.563206	...	Northwest	6.415	0.000000	27.208	0.0	24.500756	4	0	3.0	Moderate rain
2242	2015-10-16 17:00:00	Moderate rain	276.242514	3.092514	271.085195	272.543725	-0.187871	2.778891	12.095491	306.064377	...	Northwest	6.515	0.000000	27.471	0.0	25.055266	4	0	3.0	Moderate rain
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
56604	2021-12-28 04:00:00	Mild snowfall	270.414644	-2.735356	263.582322	267.519393	0.158137	2.467546	16.219425	214.752551	...	Southwest	2.007	53.587745	39.138	0.0	20.613530	5	0	5.0	Mild snowfall
56605	2021-12-28 05:00:00	Mild snowfall	270.540634	-2.609366	263.824712	267.697892	0.125990	2.431944	16.219425	211.899440	...	Southwest	1.807	53.587745	38.076	0.0	21.978703	5	0	5.0	Mild snowfall
56606	2021-12-28 06:00:00	Mild snowfall	270.638388	-2.511612	264.198556	267.826681	0.097754	2.412873	16.219425	209.410100	...	Southwest	1.509	53.587745	32.855	0.0	23.830981	5	0	5.0	Mild snowfall
56607	2021-12-28 07:00:00	Mild snowfall	270.714928	-2.435072	264.198556	267.913728	0.076540	2.405941	16.219425	206.407704	...	Southwest	1.313	53.587745	26.111	0.0	27.455241	5	0	5.0	Mild snowfall
56608	2021-12-28 08:00:00	Mild snowfall	270.801398	-2.348602	264.705794	268.001076	0.086470	2.420409	16.219425	201.798103	...	South	1.113	53.587745	17.846	0.0	26.719272	5	0	5.0	Mild snowfall

label2	wep
0.0	Moderate snowfall
1.0	Heavy snowfall with accumulated snow
2.0	Storm with freezing rain / Heavy snow- and ice pellets
3.0	Moderate rain
4.0	Snowstorm with high precipitation
5.0	Mild snowfall
6.0	Frondurchlauf / Continuous freezing rain

Draft analysis¶

Weather Data Analysis: A Regression and Classification Approach on the ERA5 Dataset¶

1. Introduction¶

1.1 Introduction to Subject - Motivation¶

1.2 Motivation for Research Question¶

1.3 General Research Question¶

1.4 Hypotheses Regarding Question¶

2. Setup¶

2.1 Methodological Insights and Impact¶

2.2 Methodology Overview¶

3. Data¶

3.1 Import data¶

3.2 Data structure¶

3.2.1 Data Description¶

3.2.2 Data Characterisitcs¶

3.2.3 Variables / Parameters¶

3.2.4 Population and sample¶

3.2.5 Overview of data and structure¶

3.2.6 Description of how the Data was collected¶

3.3 Data corrections¶

3.4 Variable lists and overview¶

3.4.1 Data Dictionary¶

3.4.2 Response Variable: Temperature¶

3.4.3 Response Variable: Weather Event¶

3.5 Data splitting¶

3.5.1 Strategies for Data Splitting¶

3.5.2 Challenges and Best Practices¶

4 Analysis¶

4.1 Descriptive statistics¶

4.2 Exploratory data analysis¶

4.2.1. Time series analysis of weather data for the region Bancroft¶

4.2.2 Distribution, Frequency and Details for Weather Parameter Observations¶

4.2.3. Further Examination of WEP and Temperature Response Variables¶

Weather Condition Frequencies¶

Overall Weather Event Distribution:¶

Extreme Weather Event Distribution:¶

4.3 Relationships¶

4.4 Plotting Data in a 3D space¶

5 Model¶

5.1 Overview Classification and Regression Algorithms¶

5.2 Potential Challenges (Feedback Project Proposal)¶

5.3. Regression Analysis:¶

5.3.1. Regression Analysis Temperature and Wind Modeling:¶

5.3.1.1 Regression Analysis Temperature and Wind Modeling: Select Model¶

5.3.1.2 Regression Analysis Temperature and Wind Modeling: Training and validation¶

5.3.1.3 Regression Analysis Temperature and Wind Modeling: Fit model¶

5.3.1.4 Regression Analysis Temperature and Wind Modeling: Evaluation on test set¶

Model Evaluation Results¶

5.3.1.6 Regression Analysis Temperature and Wind Modeling: Save Model¶

5.3.2 Linear Regression Analysis: Multiple Predictors¶

5.3.2.1 Feature Selection¶

5.3.2.2 Multiple Linear Regression¶

5.3.3 Temperature Forecast¶

5.3.3.1 SARIMAX MODEL¶

5.3.3.2 XGBOOST¶

5.3.4. Temperature Prediction (Temporal):¶

5.3.4.1 Temperature Prediction (Temporal): Select model¶

5.3.4.2 Temperature Prediction (Temporal): Training and validation¶

5.3.4.3 Temperature Prediction (Temporal): Fit model¶

5.3.4.4 Temperature Prediction (Temporal): Evaluation on test set¶

5.3.4.5 Temperature Prediction (Temporal): Visualize results¶

5.3.4.6 Temperature Prediction (Temporal): Save model¶

5.3.5 Predicting Extreme Weather Events by Temperature (Logistic Regression)¶

5.3.5.1 Predicting Extreme Weather Events by Temperature (Logistic Regression) - Select Model¶

5.3.5.2 Predicting Extreme Weather Events by Temperature (Logistic Regression) - Training and Validation¶

5.3.5.3 Predicting Extreme Weather Events by Temperature (Logistic Regression) - Fit Model¶

5.3.5.4 Predicting Extreme Weather Events by Temperature (Logistic Regression) - Evaluation on test set¶

5.3.5.5 Predicting WEP by Temperature (Logistic Regression) - Visualization¶

5.3.5.6 Predicting WEP by Temperature (Logistic Regression) - Save Model¶

5.4 Classification Analysis¶

5.4.1. Binary Classification¶

5.4.2 Multiclass Classification¶

6 Conclusions¶

Overview¶

Regression Analysis Findings¶

Classification Analysis Findings¶

Critical reflection and outlook¶